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ABSTRACT 

We study the three-dimensional (3D) hydrodynamics of the post-core-bounce phase of the collapse of a 
21-Mq star and pay special attention to the development of the standing accretion shock instability (SASI) 
and neutrino-driven convection. To this end, we perform 3D general-relativistic simulations with a 3-species 
neutrino leakage scheme with neutrino heating. Unlike "light-bulb" heating/cooling schemes, the leakage 
scheme captures the essential aspects of neutrino cooling, heating, and lepton number exchange as predicted 
by r adiation-hydrodynamics simulations. The 27-Mq progenitor was studied in 2D by B. Miiller et al. (2012; 
I arXiv: 1205.7078 ), who observed strong growth of the SASI while neutrino-driven convection was suppressed. 
In our 3D simulations, neutrino-driven convection grows from numerical perturbations imposed by our Carte- 
sian grid. It becomes the dominant instability and leads to large-scale non-oscillatory deformations of the shock 
front. These will result in strongly aspherical explosions without the need for large-scale SASI shock oscilla- 
tions. Low-^-mode SASI oscillations are present in our models, but saturate at small amplitudes that decrease 
with increasing neutrino heating and vigor of convection. Our results suggest that once neutrino -driven convec- 
tion is started, it is likely to become the dominant instability in 3D. Whether it is the primary instability after 
bounce will ultimately depend on the physical seed perturbations present in the cores of massive stars. The 
gravitational wave signal, which we extract and analyze for the first time from 3D general-relativistic models, 
will serve as an observational probe of the postbounce dynamics and, in combination with neutrinos, may allow 
us to determine the primary hydrodynamic instability. 

Subject headings: gravitation - gravitational waves - hydrodynamics - neutrinos - Stars: supernovae: general 



1. INTRODUCTION 



|Baade & Zwicky| ( 1934[ l inaugurated core-collapse super- 
nova theory with their seminal prediction that "a super-nova 
represents the transition of an ordinary star into a neutron 
star." The very basics of this theory, summarized authorita- 
tively by Bethe (1990), wer e confirmed by the observation of 
neutri nos from SN 1987 A ( |Hirata et aL]|1987| [Bionta et aL| 
|1987| l: The electron-degenerate core of a massive star (with 
mass M ^ 8 - 130Mq at zero-age main sequence [ZAMS]), 
once having reached its effective Chandrasekhar mass, be- 
comes radially unstable. Collapse ensues and, once fully dy- 
namic, separates the core into the homologous, subsonically 
contracting inner core and the outer core, which is superson- 
ically infalling. When the inner core reaches nuclear density, 
the nuclear force, which is repulsive at short distances, leads 
to a stiffening of the nuclear equation of state (EOS). The dra- 
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matically increased pressure support stabilizes the inner core, 
which overshoots its new equilibrium, then rebounds into the 
still infalling outer core. This core bounce launches a hydro- 
dynamic shock wave, which, endowed with the kinetic energy 
of the inner core, plows into the outer core. Its progression is, 
however, soon muffled by energy losses to the dissociation of 
heavy nuclei and to electron capture neutrinos that are cre- 
ated and stream out from now optically-thin regions behind 
the shock. The hydrodynamic shock thus succumbs to the ex- 
treme ram pressure of the outer core and turns into a stalled 
accretion shock. In the commonly accepted picture of the neu- 
trino mecha nism (Wilson 1985; Bethe & Wilson 1985; Bethe] 
1990[ IJanka et al.^2007^) , the shock is revived by the deposi- 
tion or neutrino energy in a layer of net neutrino heating (the 
gain layer) below the shock. In an alternative scenario, requir- 
ing very rapid progenitor rotation and efficient magnetic field 
amplification, a m agnetorotational explosion may occur (e.g., 
[Burrows et al.|20 07 and references therein). In order to leave 
behind a slowly cooling neutron star and not a black hole, 
shock revi val must occur within a few hundred millis econds 
of bounce ( [O'Connor & Ott|201 l[[Ugliano et al.]20T2 1. 



While the general picture of core-collapse supernova theory 
may be well established, details of the explosion mechanism, 
its dependence on precollapse conditions and input physics, 
and its neutrino and gravitational wave signal^remain to be 
determined by detailed first-principles numerical simulations. 



Both neutrinos and gravitational waves may be direct probes of progen- 
itor properties, supernova dynamics, and of the explosion mechanism. See, 
e.g. jOtt (2009 1; Lund et al. (2010 2012 1; Dasgupt a et al.| ( |2010| ; |Brandt etaT] 
" ^ogue et al. ( 2012j ; O'Connor & Ott ( lUl'I^ . 
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In the case of the neutrino mechanism, modern spherically- 
symmetric (ID) simulations with full Boltzmann neutrino 
transport have shown that neutrino heating alone fails to 
drive an explosion in all but the lowest-mass massive stars 
JLiebendorfer et al. 2001 , Rampp & Janka 2002; Thompson 
et al. 2003; Sumiyoshi et al. 2005, Kitaura et al. 2006, Hude- 



pohl et al.|2010 1. Spherical symmetry, however, is a poor ap- 
proximation to the situation after core bounce, even if the ini- 
tial conditions are nearly spherically symmetric. The weaken- 
ing shock leaves behind a negative entropy gradient, which is 
expected to lead to convective instability within milliseconds 
after bounce (prompt convection). Somewhat later, neutrino 
heating establishes a negative entropy gradient in the gain re- 
gion, leading to neutrino-driven convection. Strong delep- 
tonization near the neutrinosphere (where the neutrino optical 
depth ^ I; located at the edge of the protoneutron star at 
~ 10" - lO'^gcm"-') establishes a negative lepton gradient, 
driving protoneutron star convection. 



The first full axisymmetric (2D) simulations ( Herant et al.l 
T9941 [B urrows et al.|1995l[Janka & Muller|1995| |1996; Fryer 



& Heger||2000/ showed that 2D neutrino-driven convection 
could increase the efficacy of the neutrino mechanism by in- 
creasing the residence time of accreted material in the region 
of net neutrino heating and leading to high-entropy turbulent 
flow that aids shock expansion. 

A new instability, the stan ding accretion shock instability 
(SASI), was discovered by B londin et al.| ( |2003 ), who car- 
ried out idealized 2D simulations of an accretion shock using 
an analytic EOS, neutrino cooling, but no neutrino heating. 
In 2D, the SASI leads to large scale, low-order (in terms of 
spherical harmonics, £= {1,2}) deformations of the shock 
front that vary in time in a predominantly t = \ sloshing- 
type motion up and down the symmetry axis. These aspher- 
ical motions lead to larger average shock radii, increase the 
dwell time of material in the gain region, may lead to sec- 
ondary shoc ks, and are th us gene rally aiding the exp l osion 
mechanism (|Ohnishi et al. 2006; Scheck et aT]|2006| |Mur-| 
phy & Burrows 2008; Ott et al. 2008; Marek & Janka 2009', 
Miiller et al.,,2012b) . In 3D, nonaxisymmetric modes 
{m = ...,0, ...,£}) are excited as well, leading to more 
complex dynamics and smaller saturation amplitudes for in- 
dividual modes (Iwakami et al. 2008). In some 3D simula- 
tions, in particular in those that include some initial rotation, 
a strong spiral mode (£ = 1,ot = ±1), capable of redistribut- 
ing angular momentum, has been ob served ( Blondin & Mez- 
zacappa"20 07l |Iwakami et al.||2008| [2009, . Fernandez||2010i 
Wongwath anarat et al.|2010[|Rantsiou et al.j2(JTT| [ 

Perturbation theory and carefully controlled numerical ex- 
periments suggest that the SASI is driven by an advective- 
acoustic cycle in which entropy and vorticity perturbations are 
advected from the shock front to the edge of the protoneutron 
star. There they trigger the emission of acoustic perturbations 
that travel upstream in the subsonic flow of the postshock re- 
gion and amplify perturbations in the shock front, thus creat- 
ing a feedback cy cle that injects power preferentially into low- 
order modes (s ee Foglizzo 2002 ; Foglizzo et al. 2006, 2007"; 
Ohnishi et al.ll2006; Yamasaki & Yamada 2007; Fernandez, 



& Thompson 2009b a; Scheck et al. 2008 , Guilet & Foglizz^ 
5Dl2|and references therein). The saturation of the SASI has 



In a real core-collapse supernova, neutrino-driven convec- 
tion and SASI overlap in space and may grow at the same 
time. In the linear regime, in which seed perturbations are 
minute, they can be clearly separated: SASI's fastest growing 
mode is 1= 1, while convective eddies will grow with hori- 
zontal wavelengths a few times of the entropy scale height, 
giving ^ ~ 7 - 8 in the postbounce supernova context ( Herant 



|et al. 1992[ |Chandrasekhar|[T96 1 | l . In convection, however 



all modes are unstable and will eventually grow to nonlinear 
amplitudes if convection is able to develop at all. 

How convection and SASI interact in the nonlinear regime, 
which of them becomes the dominant instability, how this 
may depend on the dimensionality (2D vs. 3D), and the ram- 
ifications of all this for the explosion mechanism are open 
questions that are currently under much debate. 



been proposed to occur via parasitic Rayleigh-Taylor and/or 
Kelvin-Helmholtz instabilities that operate on the entropy gra - 
dients and vorticity generated by the SASI ( [Guilet etal.|2010| l. 



Foglizzo et al. 1(2006) argued, based on linear theory, that in 
the absence of large (i.e., nonlinear) perturbations the devel- 
opment of neutrino-driven convection may be suppressed if 
slowly developing eddies are advected out of the convectively 
unstable region before they can grow significantly. In this sce- 
nario, SASI would be the primary instability . This was also 
found in the idealized simulations of Oh nishi et al.| ( |2006[ ), 
who studied the 2D evolution of an artificially set up accre- 
tion shock with a constant accretion rate and analy tic neu- 
trino cooling , heating, and deleptonization functions. Scheck 
et al. (2008) performed 2D energy-averaged (gray) neutrino 
radiation-hydrodynamics postbounce simulations of a \5-Mq 
progenitor star in a carefully controlled setting to study the 
develo pment of the SA SI. They, too, confirmed the result 
of Foglizzo et al. ( 2006| l and showed that if sufficiently large 
(> 1%) perturbations from sphericity are present in the up- 
stream flow, neutrino-driven convection becomes the primary 
and dominant instability. 

If linearly-growing convection is suppressed by high ad- 
vection velocities in the gain region, then one would ex- 
pect a dependence of the relative importance of SASI and 
convection on the postbounce accretion rate and, hence, on 
the progenitor s tar. This was convincingly confirmed by the 
recent work of |B. Miiller et al.[ ( [2012a ), who carried out 
full first-principles 2D general-relativistic (GR) multi-energy 
radiation-hydrodynamics postbounce simulations of a 8.1- 
Mq low-metallicity star with a small core and low postbounce 
accretion rate and of a 21 -Mq star of solar metallicity with a 
large core and high accret i on rate . In agreement with the pre- 
diction of [FogUzzo et al.[ ( [2006| l, they found strong convec- 
tion and absent SASI in the 8.I-M0 star and strong SASI and 
nearly absent convection in the 21-Mq progenitor. In both 
cases, explosions developed within ^200 ms of bounce. 

In a different line of research targeted at understanding the 
dependence of the neutrino mechanism on dimensionality, 
[Nordhaus eraL] ('2010) carried out ID, 2D, and 3D Newto- 
nian collapse simulations of a 15-Mq progenitor They used 
the simple analytic heating and coo ling prescription intro- 
duced by Murphy & Burrows[ ( p008) ) (hereaf ter the MB08 
"light-bulb" scheme) on the basis of the work of Janka ( 2001 1. 
Their 3D simulations did not show a dominant £ = 1 oscilla- 
tory SASI m ode observed in 2D (Schec k et al.|2006[phnishi| 
[et al. 2006; [Murphy & Burrow s 2008). Using the critical 
lumino sity vs. accretion rate approach of Burrows & Goshy 
( 1993 1, they reported that in 3D explosions could be obtained 
at ^-^15-25% and ~'40-50% lower neutrino luminosities than 
in 2D and ID, respectively. 
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The |Nordhaus et aL](|2010| l 3D vs. 2D result was not con- 
firmed by Hanke et al. ( 2012[ l. These authors performed New- 
tonian simulations with neutrino approximations very similar 
to the MB08 light bulb, but used a different 3D hydrodynam- 
ics code. They did not find clear evidence that 3D effects 
facilitate the development of an explosion to a greater degree 
than the non-radial motions due to SASI and convection in 



2D. However, in agreement with |Nordhaus et al. (2010), they 
did not find large-scale oscillatory low-order modes in their 
3D simulations. They hypothesized that this may be less of 
a 3D effect than an effect of the rather simple treatment of 
neutrino heating and cooling by Nordhaus et al. (2010). The 
arguably greatest limitation of the MB08 light-bulb scheme is 
its inability to track the contraction of the protoneutron star, 
leading to too low advection velocities in the gain region, thus 
artificially favoring neutrino-dr i ven co nvection over the SASI. 
The results of Takiw aki et al.| ( |2012[ ), whose Newtonian 3D 
simulations used a multi-energy approximate neutrino trans- 
port scheme, appear supportive of this assertion. However, 
these simulations were carried out with very low resolution 
and the low-order modes appear to be clearly oscillatory only 
at early times. 

Using the same MB08 light-bulb approximation for neutri- 
nos and an updated version of the Nordhaus et al. (2010) code, 
IBurrows et al. (2012), Murphy et al. (2012 ), and Dolen ce etal.| 
(2012 1 performed and analyzed another set of 2D and 3D sim- 
ulations to investigate the roles of SASI and convection in 
the postbounce evolution of a 15-Mq progenitor. Compar- 
ing 2D and 3D resul ts for the ev olution of low-order fluid 
mode amplitudes. Burrows et al. (2012) showed that at the 
same MB08 driving luminosity, oscillatory mode amplitudes 
are much smaller in 3D than in 2D. In models that develop 
an explosion a non-oscillatory £ = I dipole asphericity grows 
already in the early postbounce evolution. Furthermore, they 
showed that the oscillatory £=l modes observed in 2D - and 
generally associated with the SASI - occur even in the case of 
a high light-bulb driving luminosity, in which neutrino-driven 
convection is the dominant instability. They argued that in 
successful explosions by the neutrino mechanism, neutrino- 
driven convection should be the dominant instability . How- 
ever, for the reasons put forth by |Hanke et al.| ( |2012[ ) and [BT] 



IMiiUer et al.| ( |2012a i and discussed in the above, the predic 
five power or these Iight-bulb simulations may be limited. 

Ultimately, high-resolution 3D energy-dependent GR neu- 
trino radiation-hydrodynamics simulations will be needed for 
final answers regarding the explosion mechanisms and the 
role of the various instabilities involved. Such simulations are 
computationally extremely challenging and current attempts 
are forced to use low spatial resolution (Takiwaki et al.|2012( 



Kuroda et al. 2012 ), the gray approximation (Wongwathanarat 
et al.|2010,,E. Muller et al..2012,,Kurod a et al. 2012;, and/or 



employ an arti ficial inner boundary, cu tting out t he protoneu 
tron star core (IWongwathanarat et al.||2010£ lE. Muller et al. 

m^ - — 

In this paper, we present results from 3D hydrodynamic 
postbounce supernova calculations that attempt to strike a bal- 
ance between the computationally cheap, but possibly too 
simplistic light-bulb approximation and true 3D radiation- 
hydrodynamics simulations, which cannot yet be performed 
without at least partially debilitating limitations. Our simu- 
lations use the Zelmani core collapse simulation package 
( |Ott et al^|2012| i and are fully general relativistic. We make 



no symmetry assumptions and use no artificial inner bound- 
ary. We employ a novel computational setup with a multi- 
block approach that provides curvilinear grid blocks to track 
the collapse of the outer core and Cartesian adaptive-mesh re- 
finement (AMR) grids covering the central region, including 
the protoneutron star and the entire shock. We treat neutri- 
nos in the postbounce phase with an energy-averaged three- 
species neutrino leakage scheme with neutrino heating. The 
only free parameter of this scheme is a scaling factor in the 
charged-current energy deposition rate. As we shall demon- 
strate, the leakage scheme captures the essential aspects of 
neutrino cooling, neutrino heating, and lepton number ex- 
change and yields a much more realistic postbounce configu- 
ration than possible with the MB08 light-bulb approach. 

We apply Zelmani to the collapse and postbounce evolu- 
tion of the 27-Mp) pr ogenitor star that was considered by'BT] 



Muller et al. (2012a I and shown to be highly susceptible to 
the SASI in their fully self-consistent 2D GR simulations. B. 



Muller et aL] ( j2012a| find a SASI-aided explosion that deveT 
ops within ~ 150 -200 ms after bounce, making this progeni 
tor ideal for studying the SASI in computationally expensive 
high-resolution 3D simulations. We carry out four simula- 
tions of the 21-Mq progenitor, varying the strength of neu- 
trino heating. We evolve these four models from the onset of 
collapse to ^150- 190ms after bounce at an effective angular 
resolution of 0.85° at a radius of 100km. The linear resolution 
at this radius is ^1.5 km. The maximum resolution covering 
the protoneutron star core is ^370 m. 

We find that neutrino-driven convection is able to grow 
from the numerical seed perturbations imposed by our Carte- 
sian AMR approach. It becomes the dominant instability in 
the postbounce dynamics of all of our models. In the case of 
strong neutrino heating, convection, which is initially man- 
ifest as small-scale cells of rising hotter and sinking cooler 
material, develops into large blobs of high entropy mate- 
rial. These push out the shock and lead to large-scale non- 
oscillatory shock deformations. We also observe growth of 
oscillatory low-(^,m) deformations associated with the SASI. 
However, these saturate at small amplitudes that decrease fur- 
ther with increasing strength of neutrino heating and vigor of 
convection. The SASI remains sub-dominant at all times in 
our simulations. Our results suggest that if neutrino-driven 
convection is able to grow in 3D - which will generally de- 
pend on the postbounce ac cretion rate an d on the seed pertur- 
bations present in the flow (jScheck et al. 2008 ) - it will dom- 
inate the postbounce flow. We extract the gravitational wave 
signals generated by accelerated quadrupole mass motions in 
our models and find that the strongest component of the signal 
comes from the initial burst of convection, which grows on the 
negative entropy gradient left behind by the stalling shock. 

This paper is structured as follows. In Section |2] we de- 
scribe Zelmani and give details on grid setup, EOS, the 
leakage/heating scheme, and the progenitor model. In Sec 
tion[3]we present the results of our simulations. First, in 5 3. 1 



we give an overview of the overall postbounce evolution of 
our models. We then discuss in detail the postbounc e co nfig- 

the 



urations resulting from our leakage/heating scheme (9372 



development of neutrino-driven convection and SAS I {\ 3.3 1, 
various criteria for neutrino-driven explosions ( ^3.4[ i, an d th e 
gravitational wave signals extracted from our models ( ^3.5[ ). 
We summarize our findings and conclude in Section]?] 
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2. METHODS AND INITIAL CONDITIONS 

We carry out our 3D GR simulations with the Zelmani 
core collapse simulation package. Zel mani is based on th e 
open-source Einstein Toolki tp](|Loffler et al.||2012|l 



for numerical relativity and relativistic computation al astro- 
physics. It builds upon the Carpet AMR driver (Schnet- 



ter et al.| 


2004; 


ney et al. 


120111 



J 



within the Cactus 



2.1. Spacetime Evolution and Hydrodynamics 

We evolve the full Einstein equations without approxima- 
tions in a 3 + 1 decomposition as a Cauchy initial boundary 
value problem (see, e.g., Baumgarte & Shapiro 2010), us- 
ing the conformal-traceles s BSSN formulation (Baumgarte & 
Shapiro] 19991 IShibata & Nakamur a^l995) . A 1 + log slicing 



condition ( [Alcubierre et al. 2000) controls the evolutio n of 



the lapse function a , and a modified F-driver condition ( |A1- 
[cubierre et al. 12003] ) is used for the evolution of the coordinate 
shift vector /3 . The BSSN equations and the gauge conditions 
are implemented in the module CTGamma using fourth-order 
ac curate finite d i fferen cing. Implementation detai ls are given 
Pollney et aLl ( |2011| l and |Reisswig et al.| ( |20T2] i. 



m 



We use a flux-conservative formulation of the GR Eu- 
ler equations, implemented in the GR hydrodynamics mod- 
ule GRHydro, which is part of the Einstein Toolkit. 
GRHydro is based on a finite-volume high-resolution shock- 
capturing scheme and works with general finite-temperature 
microphysical EOS. We employ the enhanced piecewise- 
parabolic method for reconstruction of state variables at cell 
interfaces ( [McCorquodale & Colella|[2011[ [Reisswig et al.| 
2012| i and subsequently solve approximate Riemann prob 



lems to compute intercell fluxe s with the HLLE solve r ( Ein- 
[feldt 1988 ). Details are given in Reisswig et al. (20121. 



Both spacetime evolution and GR hydrodynamics are dis- 
cretized in a semi-discrete fashion and coupled with the 
Method of Lines (Hyman 1976 ) using a multi-rate Runge- 
Kutta integrator (Reisswig et al. 2012 ), providing fourth-order 
and second-order accuracy in time for spacetime and GR hy- 
drodynamics, respectively. The time step is limited by the 
speed of light and we use a constant Courant-Friedrichs-Lewy 
factor of 0.4. 

2.2. Multi-Block Infrastructure, Adaptive Mesh Refinement, 
and Grid Setup 



We employ th e multi-block infr astructure Llama ( Pollney 
|et al.p Oll , Rei sswig et al.||20T2| , which allows us to cover 
the computational domain using a set of overlapping curvi- 
linear grid blocks that are logically Cartesian but physically 
curvilinear (so-called "inflated cubes"), adapted to the overall 
spherical topology of the collapse problem. We employ a set 
of such curvilinear blocks to track the collapse of the outer 
core, while the interior domain containing the protoneutron 
star, the postshock region, and the shock itself is covered by 
an adaptively refined Cartesian mesh. 

The spherical inflated-cube multi-block system discretizes 
one spherical shell via six angular grid blocks designed 



http : //www . einsteintoolkit . org 



such that one angular coordinate direction always coincides 
at inter-block boundaries. This allows us to use efficient 
fourth-order one-dimensional interpolation to update ghost 
zone information between neighboring blocks ( Thornburg] 
[2004 | l. Furthermore, this particular multi-block system of- 
fers an almost uniform distribution of points across the 
sphere (i.e. without clustering of points at the poles), thus 
avoiding distortions and pathologies associated with standard 
spherical-polar grids. 

The adaptively refined central Cartesian block is based 
on cell-centered and flux-conservative mesh-refinement tech- 
niques, provided by the open-source AMR driver Carpet 
(Schn etter et al.l 2004 ; Reisswig et al |2012| l. AMR is imple- 
mented with subcycling in time, following the approach of 
Berger & Oliger (1984). We make use of refluxing, which 
correctly adjusts fluxes at mesh refinement boundarie s after 
the AMR restriction operation (Reisswig et al. 2012). This 
ensures that mass, momentum, and energy fluxes are ex- 
actly conserved, even in the presence of strong shocks and 
other discontinuities. To update zones at AMR boundaries 
and to initialize new grid points after regridding, we make 
use of fourth-order prolongation for the spacetime curvature 
variables, and second-order essentially non-oscill atory pro- 
longation for the matter variables. As detailed in [Reisswig] 
|et al.| (j2012 ), spacetime variables are restricted from fine onto 
coarse grids using a third-order polynomial, while matter vari- 
ables are restricted via cell averaging. 

All simulations are carried out with the same general grid 
setup. In the central region, we use five nested Cartesian 
grids with a factor of 2 in resolution between each of them. 
The finest grid has a linear cell size dx = 0.37km and ex- 
tends out to 17.7km. The second finest grid has a cell size 
of dx = 0.74km and extends to 59 km, while the third grid has 
dx= 1 .48 km and is set up to adaptively track the shock, en- 
suring that shock itself and the turbulent flow in the gain layer 
behind the shock are always resolved with no worse resolution 
than dx = 1.48 km. For a shock radius of 100 km, this cor- 
responds to an effective angular resolution dx/R of ^0.85°. 
There are two additional coarser grids with dx = 2.95 km and 
dx = 5.9km in the Cartesian region, which extends to 532km, 
where it overlaps with the outer spherical cube grid. The lat- 
ter's radial cell size dr at its inner boundary is the same as 
the dx of the coarsest Cartesian grid it overlaps with, dr is 
held constant out to a radius of ^3000 km and then smoothly 
reduced to dr = 189km at the outer boundary at ~'15000km. 
Each of the six cubed-sphere blocks has 3 1 angular zones each 
in polar angle 9 and azimuthal angle ip. This corresponds to 
an effective cell size of ^2.9°. 

We start our simulations at the onset of collapse with only 
the coarsest of the Cartesian AMR grids active and progres- 
sively activate the finer grids when the central density in the 
collapsing core reaches 3.2 x 10" gem"-', 1.3 x lO'^gcm"-', 
5.1 X lO'^gcm"-', and 2.0 x lO'-'gcm"-', respectively. 

2.3. Equation of State 

We employ a tabulated versio n of the finite-temperature nu- 
clear EOS by |Lattimer & Swes ty ( 1991 ). This EOS is based 
on the compressible liquid-drop model with a nuclear sym- 
metry energy of 29.3 MeV. We use its variant with a nuclear 
compression modulus Kq of 220 MeV, since it yields a cold 
neutron star mass-radius relationship in agreement with cur- 
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rent observational and theoretical constraints (e.g., Demorest 
let al.|2010[|Hebeler et al.|20T0l|Steiner et al.|2010| i. 

We employ the Lattimer-Swesty EOS at densities above 
lO^gcm"^, where T > 0.5MeV at all times in the core col- 
lapse context and nuclear statistical equilibrium ( NSE) holds. 
At low er densities, we employ the Timmes EOS (Timmes & 
Arnett|199 9 ) and assume that the matter is an ideal gas com- 
posed of electrons, positrons, photons, neutrons, protons, al- 
pha particles, and heavy nuclei with the average A and Z given 
by the Lattimer-Swesty EOS at the transition density. This 
is an approximation and may lead to slightly incorrect pres- 
sures in non-NSE regions that result in changes in the col- 
lapse times for the silicon and carbon/oxygen shells. Ideally, 
a fully consistent treatment with multiple advected chemical 
species, a nuclear reaction network and tr ansition in and out of 
NSE with a NSE network as proposed by Buras et aL] ( |2006b| 
should be implemented. This, however, is beyond the scope 
of the present study. 

Details on the EOS table and on the implementation of the 
contribution of electrons, positrons, and photons, as well as 
ot her details of the construction of the table are described 
in [O'Connor & Ott| ( |20T0 ). The table itself as well as table 
generation and interpolation routines are available at |http : | 
|/ / www .stellarcollapse. org, 

2.4. Neutrino Treatment 

We employ the approximate neutrino treatmen t of the 
open-source code GRID (O'Co nnor & Ott] 20101, which 
was adapted to 3D and imple mented m the module 
ZelmaniLeak by|Ott et al.|(|2012|l. The source code is avail- 
able from jhttp : //www, stella rcolla pse . org 

Before core bounce, the primary neutrino emission process 
is electron capture on free and bound protons, leading to a 
reduction of the electron fraction in the collapsing core. 
We include this effect and associated changes of the specific 
entropy in the approximate way proposed by Liebendorfer] 
( [2005 p . He showed, on the basis of ID Boltzmann neutrino 
radiation-hydrodynamics simulations, that Yg in the collapse 
phase can be well parameterized as a function of rest-mass 
density p. This parameterization shows only small variations 
with progenitor star and nuclear EOS. We employ an ana- 
lytic Yeip) fit to the results of ID radiation-hydrodynamics 
collap se simulations of a 20- Mq solar-metallicity progenitor 
star of'Woosley et al.'("2002^ obtained with the code and mi- 
crophysics of Buras et al. ( 200 6b|l. The same Ye{ p) profile was 
used in |Ott et al.| ( |2007a|bj i and [Ott et aL] ( [2012| i. 

In the late collapse phase, when neutrinos begin to be 
trapped in the inner core, and throughout the postbounce 
phase, momentum exchange between neutrinos and matter be- 
comes non-negligible. The effect of this "neutrino stress" is 
naturally captured by the coupling of radiati on and matter in 
full n eutrino transport calculations (see, e.g., B. Miiller et al. 
[2010| l. In our approximate treatment, we must include it ex- 
plicitly. We assume that neutrino stress is relevant only above 
a fiducial trapping density of 2 x lO'^gcm"-' and approximate 
the stress as the gradient of the neutrino Fermi pressure. The 
stress is then included as a source term in the GR hydro- 
dynamics equations at each time-integration substep and the 
neu trino Fe rmi pressure is included in the stress- energy tensor 
(see jOtt et al. 2007b, and .O'Connor & Ott|2010| for details). 



After core bounce, which we define as the time at which 
the specific entropy at the edge of the inner core reaches 
Sfefibaryon"', signaling shock formation, the simple Y^p) ap- 
proximation breaks down and fails to even qualitatively cap- 
ture the effects of neutrino processes occurring in the post- 
bounce phase. Dissociation of iron-group nuclei by the shock 
provides a sea of free protons for electrons to capture on, lead- 
ing to the neutronization burst of electron neutrinos (ve) and 
a steep drop of Ye in the region just outside the nascent pro- 
toneutron star. High temperatures and low Ye in the lower 
postshock region allow for the appearance of positrons that 
capture on neutrons, leading to the emission of electron an- 
tineutrinos {Ve)- High temperatures in the protoneutron star 
core lead to neutral-current pair emission of neutrinos of all 
species. 

In order to capture the aforementioned processes and their 
effects in terms of cooling, heating, and deleptonization in the 
region behi nd the shock, we switch to the neutrino leakage 
scheme of O'Connor & Ott (2010) at bounce. We consider 
three neutrino species, Vg, Vg, and v^ = {Vf^^v^^v-nVr}, where 
we lump the heavy-lepton neutrinos together, since they par- 
ticipate only in neutral current processes and have very similar 
cross sections in the core-collapse supernova environment. 

The leakage scheme provides approximate energy and 
number emission and absorption rates based on local thermo- 
dynamics and the optical depth in the postshock region. Neu- 
trino absorption and emission are ignored outside the shock. 
The optical depth requires a non-local calculation, which we 
solve in a ray-by-ray way, computing an optical depth integral 
r,^, along radial ra ys cast into 9 and directions (see Fig. 1 
of Ott et al.]|2012 1 from the origin. We then interpolate tri- 
linearly in (r,0,(p) to obtain the optical depth at the centers 
of Cartesian grid cells. Ideally, an optical depth calculation 
should be carried out into all directions from any given cell 
and the minimum v alue should be used as the optical depth 
of that cell (see, e.g., |Ruffert et al.|1996p . However, for situa- 
tions that are spherical at zeroth order, like the one considered 
here, the computationally much cheaper ray-by-ray approach 
should be sufficient. In our simulations, we employ 37 rays 
in 9, covering [0,7r], and 75 rays in (p, covering [0,27r]. Each 
ray has 800 equidistant points to ^ 600km and 200 logarith- 
mically spaced points covering ^ 600 -3000 km. 

We approximately include neutrino heating by charged- 
current absorption of h'e and on neutrons and protons, re- 
spectively. For this, we ma ke use of a lo cal heating function 
based on the derivations by Janka (2001 1, 



-2r„ 



(1) 



Here L^. is the neutrino luminosity incident from below, 5^ = 
0.25(1 + 3Q!^)(Jo(mgC^)"^, where (Tq is the fiducial weak inter- 
action cross-section ^ 1.76 x 10~'''*cm^, a = 1.23, and nieC^ is 
the electron rest mass energy in MeV. p is the rest-mass den- 
sity, m„ is the neutron mass in grams, and Xj is the neutron 
(or proton) mass fraction, (ej, ) is the mean-squared energy 
of Vi neutrinos. We approximate it by taking the matter tem- 
perature TiMS iy, at the neutrinosphere (where t^. = 2/3) and 
evaluating 



= 7; 



NS, 



-^5(?yt/„Ns) 



(2) 



where the J^j are Fermi integrals J-j(r]) = /„°° dxx^e^ + 1) 
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Table 1 

Key Simulation Parameters and Results 



Model 


/heat 




de,d(j) 


^end 


^shock.max 


^shock.av 


^shock.min 






(km) 


@ 100 km 


(ms) 


@'e„d 


@'e„d 


@'end 








(degrees) 




(km) 


(km) 


(km) 


i27/heatl.00 


1.00 


1.48 


0.85 


184 


82 


71 


62 


i27/hea,1.05 


1.05 


1.48 


0.85 


192 


259 


189 


152 


i27/hea,1.10 


1.10 


1.48 


0.85 


165 


428 


306 


204 


i27/hea,1.15 


1.15 


1.48 


0.85 


154 


432 


336 


267 



Note. — /iieat is the scaling factor in the neutrino heating rate (Eq.fTJ, rf-Vshock is the mini- 
mum linear resolution covering the shock and the region interior to it. aa,d(j> @ 100 km is the 
effective angular resolution at a radius of 100 km, (^nj is the time after core bounce at which the 
simulation is stopped, and Rshock.maxi '^shock.av and Rshock.min are the final maximum, average, 
and minimum shock radius, respectively. 



and i]y.^Ns = Mi^i.Ns(fcB7NS,iy,) where ^^^ns is the chemical 
potential of neutrino species Vj at its neutrino sphere. The 
factor (F~^ ) is the mean inverse flux factor, which depends on 
details of the neutrino radiation field. We parameterize it as a 
function of optical depth t^. based on the angle-depe ndent ra- 
diation fields of the neutrino transport calculations of |Ott et aT] 
^008) and set (F"') = 4.275r^, + 1.15. While the true mean 
mverse flux factor will asymptote to 1 at infinity, this simple 
fit leads to values in the postsh ock region (we in clude heating 
only there) in agreement with |Ott et al.| ( p008] l. Finally, the 
factor e"^^"/ is applied to strongly suppress heating at opti- 
cal depth above unity. The leakage scheme implementation 
in Zel mani varies sligh tly from the orig inal implementa- 
tion in lOTonnor & Ott| (2010, 20TT]l. In [O'Connor & Ott| 



(2010 201 \}, leakage was calculated only inside the shock to 



avoid unnecessary calculations outside of the shock where lit- 
tle cooling or heating occurred. To facilitate easy implemen- 
tation in Zelmani, where the angle-dependent shock radius 
is evaluated only infrequently, we have removed the explicit 
dependence on the shock radius and have replaced it with a 
condition on the mass fraction of heavy nuclei: we only cal- 
culate the heating and cooling terms where the heavy nuclei 
mass fraction is smaller than 0.5 or the density is higher than 
gcm" 



We obtain the neutrinosphere locations and the thermody- 
namic conditions for Eq. (j2]i from the rays used for the optical 
depth calculations. We also solve full leakage problems in- 
cluding heating along the rays to obtain an estimate for the 
incident luminosity L^, needed by Eq. ([T]i. The estimates for 
L^. and (e^ ) are then interpolated between rays for the local 
leakage calculations in Cartesian grid cells. 

All of the above leakage calculations are carried out 
operator-split after the fully coupled spacetime/hydro update 
and are first order in time. We find this to be sufficiently ac- 
curate and stable, due to the small time step imposed by the 
light travel time through the smallest cell. The energy and lep- 
ton number updates are applied to the fluid rest-frame quan- 
tities and we ignore velocity dependence or other relativistic 
effects in consideration of the overall very approximate nature 
of the leakage scheme. The computationally most expensive 
aspect of the leakage scheme is the interpolation of density, 
temperature, and electron fraction onto the rays. This inter- 
polation is executed at every time step in the highly dynamic 
early postbounce phase. We later switch to carrying out this 
interpolation only every 16 fine grid time steps (correspond- 
ing to every ^ 8 x 10"^ s) while continuing to evaluate the 
local expressions at every time step. 



2.5. Initial Model 

We simulate core collapse and postbounce evolution in the 
nonrotating sin gle-star 21-Mq solar-metallicity model 527 of 
[Woosley et al.| ( 2002) . We choose this particular model to fa- 
cilitate compa risons with the rece nt 2D results of B. Miill^ 
eTaL] ( |2012a i As pointed out by |B. MuUer eta\.^ ( |2012a| 7 



1 .5 Mq and a 



this progenitor has an iron-core rnas^"or 
silicon-shell mass of ^0. 18M0. AccOTding to O'Connor & 



,Ott|(201 1) , this progenitor, having a bounce compactness pa- 
rameter ^2.5 = 2.5(R[M = 2.5M0)]/lOOOkm)-' = 0.233, is a 
likely candidate for expl osion via t he neu trino mechanism. 
The recent work of Ugliano et al.| ( |2012[ ) predicts a higher 
failed core-coll apse supernova ra te for the solar metalhcity 
model set of Woosley et al. ( 2002 1 than the work of |0'Connor 



& Ott (201 Ijl. However, they also predict that this particular 
presupernova model is a progenitor of a successful neutrino- 
driven core-collapse supernova. 



Using the sphe rically-symmetric GRID code of [O'Connor 
& Ott (2010 [201 1) ) and modifying its leakage scheme to be 
identical to what we use in Zelmani, we find that fheai =1.18 



heat " 

is the critical value of the scaling factor in Eq. (T]) to drive 
an explosion that sets in at late times after multiple cycles of 
radial shock oscillations, /heat = 1-28 is required to drive an 
explosion without shock oscillations that sets in at ^ 150ms 
after bounce. 

We map model i27 to our 3D grid under the assumption 
that the ID profile data represent cell averages and use the 
radii of cell centers for interpolation. The initial spacetime is 
set up under the assumption of spherical symmetry and weak 
gravity, using the Newtonian line element without distinction 
between areal and isotropic radius. 

3. RESULTS 

3.1. Overall Postbounce Evolution 

We simulate core collapse and bounce of the 21-Mq pro- 
genitor in full 3D with adaptive mesh refinement, adding re- 
finement levels as the collapse towards a protoneutron star 
proceeds (see 5 2.2 1. Core bounce, defined as the time when 
the entropy at the edge of the inner core reaches 3 ke baryon"' , 
occurs at ^299 ms. At bo unce, we sw itch from the Yeip) pa- 
rameterization of Lieb endorfer[ p005] l to the leakage/heating 
scheme described in ^2.4[ This scheme includes a scaling 

'^^ We define the iron-core mass as the mass coordinate that has a Ye of 
0.495 
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Figure 1. Global evolution of the shock in all models. Top panel: av- 
erage shock radii (i?shock). Center panel: standard deviation (Tshock = 

[(47r)"' / iin[/?shock-(Sshock}]^] '''^ of the shock radii. Bottom panel: ra- 
tio of maximum to minimum shock radius. The shock radii of models with 
/heta ^ 1-05 exhibit positive trends in their average shock radii and have 
growing o"shock ™d ratios between maximum and minimum shock radius. 
Model i27/heatl.00's shock radius starts decreasing at ~100itis after bounce 
and its o"si,o^;, and min/max shock radii ratios oscillate around moderate val- 
ues. 

factor /heat in the neutrino energy deposition rate (Eq. [T]l. 
We carry out four long-term postbounce simulations, choos- 
ing /heat = {1.00,1.05,1.10,1.15} to study the influence of 
changes in the heating rate on the postbounce evolution. All 
models are labeled according to their value of /heat- For exam- 
ple, i27/heatl -00 is the model with /heat = 1 -00. All models are 
evolved to > 150 ms after bounce and for as long as our com- 
puter time allocations allow at a cost of ^ 25,000 CPU hours 
per millisecond of physical postbounce time (see Table[T]l. 

In the top panel of Fig. [1] we show the angle-averaged 
shock radius as a function oftime in the four simulated mod- 
els. After the early dynamic expansion phase, shock expan- 
sion stagnates and the shock stalls at 100- 130 km about 40 ms 
after bounce. Up to this point, the evolution is virtually inde- 
pendent of /heat- In the subsequent period of quasi-stationary 
evolution, the gain layer develops and neutrino heating drives 
a secular shock expansion. 

The neutrino luminosity emitted from the protoneutron star 
core and provided by accretion is identical in all models. 
Hence, as shown in the top panels of Fig. [2] there is a mono- 
tonic increase with /heat in the net neutrino heating rate Qnet 
and in the heating efficiency r; = QmiiLi,^ +Li>^)~', where we 
use the angle-averaged luminosities at the base of the gain 
layer. Varying /heat by a moderate 15% from 1.00 to 1.15 re- 
sults in ~100% more total net heating, since the increase in 
the local energy deposition rate results in an expanded gain 
layer with more mass that is able to absorb net neutrino en- 
ergy (cf. bottom panel of Fig.|2|i. 




' - 'bounce [ms] 



Figure 2. Evolution of key integral quantities indicative for the strength of 
neutrino heating. Top panel: net neutrino heating rate Q""^' (total heating 
iTunus total cooling). Center panel: heating efficiency r; defined as the net 
heating rate divided by the sum of the Ue and Ue angle-averaged luminosities 
incident below the gain layer Bottom panel: Mass in the gain layer (left or- 
dinate) and density-weighted average specific entropy in the gain layer (right 
ordinate). Heating rate, efficiency, and mass in the gain layer all increase 
inonotonically with increasing heating scaling factor /heat. Interestingly, the 
specific entropy average (igain) in the gain layer does not exhibit such a de- 
pendence on /iieat and the (.Vjai„) curves of all models are nearly identical until 
> 100 ms after bounce, at which point the overall hydrodynamic evolutions 
have diverged. 

The quantitative differences in neutrino energy deposition 
translate to qualitative differences in the shock evolution. In 
model .s27/heatl-00, which has the least heating, shock stag- 
nation turns into recession and the average shock radius de- 
creases to ^70km at the end of the simulation. The situa- 
tion is very different in models 527/heatl-lO and .s27/heatl-15, 
which both show expanding average shock radii, surpassing 
300 km at the end of their simulations and trending towards 
explosion. Model .s27/heatl-05 is somewhere in between, but 
has a slowly, but steadily increasing average shock radius that 
reaches ^190 km at the end of the simulation. 

The center and bottom panels of Fig.[T]display simple mea- 
sures of the asphericity of the shock: CTshock, the angular stan- 
dard deviation of the shock radius, and /^shock.max/^shock.min, 
the ratio of maximum to minimum shock radius. Both quan- 
tities show an initial local maximum at ^8ms after bounce, 
which is due to an initial transient large £ = 4 deformation of 
the shock front caused by the Cartesian grid employed in our 
simulations. We will discuss this further in ^3.31 In the first 
40 ms after bounce, all models show very similar small de- 
viations of the shock from spherical symmetry. Differences 
between models begin to be apparent at the same time their 
average shock radii begin to diverge. Models 527/heatl-lO 
and 527/heatl-15 exhibit very large asymmetries with CTshock ^ 
30km and almost a factor of two in radius between maximum 
and minimum shock radius at the end of their simulations 





Figure 3. 3D Volume renderings of the specific entropy at ~150ms after bounce in tlie four simulated models. The z-axis of the frames is the vertical, x is 
the horizontal and y is into the frame. The scale of the frames is 700 km on a side. The colormap is chosen such that cyan corresponds to a moderate specific 
entropy of ~4.3/:b baryon"' , indicating the shock front and low-entropy regions near the protoneutron star. Regions in yellow indicate higher entropy gas at 
.V ~ 16<:b baryon"' and red regions correspond to gas with s ~ 20A:b baryon"' . These values are chosen to highlight the surface of the shock and gas at a 
representative "intermediate" and a representative "high" specific entropy. Note the large scale global asymmetries and the many small blob-like protrusions in 
the shock fronts of models whose shock has reached large radii. 



(see TablefTlfor final minimum, maximum, and average shock 
radii for alTmodels). Model ,s'27/heatl 05 also shows growing 
asymmetry with increasing postbounce time, similar to the 
models with /heat =1-10 and 1.15, but, at least in dshock, a peri- 
odicity is visible, which is lacking completely or is occurring 
at a much smaller level in the two models with larger /heat- In 
model .s27/heatl-00, which does not show a positive trend in 
its shock radius, the deviations of the shock from sphericity 
remain small and maximum and minimum shock radius dif- 
fer, on average, by ^20% and this average difference does not 
grow until the end of the simulation. There is, however, clear 
oscillatory behavior (with a short period of ^15ms) in this 



model's shock radius variations, which may be indicative of 
SASI activity. We shall investigate this further in 5 3.3 

In 

tropy 



Fig. [3] we present volume renderings of the specific en- 
' at ~T50ms after bounce for all four models. The render- 
ings are all plotted at the same scale to emphasize the differ- 
ences in shock radius and 3D geometry between the models. 
The color map and rendering opacity are chosen to emphasize 
(0 regions with specific entropy of ~'4.3^Bbaryon~' (cyan), 
(//) regions with a representative "intermediate" specific en- 
tropy of ^16fcBbaryon~' (yellow), and, (///) regions with 
a representative "high" specific entropy of ~20^Bbaryon~' 
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527/heatl-OO 



s27/heatl-05 



527/heatl-15 




Specific Entropy 
(AiB/baryon) 



i — ^bounce — 154 ms 



Figure 4. Colormaps of the specific entropy in tlie x-z. plane in models i'27/i,e;„ 1 .00 (left column), sllfi^^.^t 1 .05 (center column), and sllfi^^^f 1.15 (right column) 
at 80, 115, and 154ms after core bounce. The linear scales of the three vertical panels are 350, 450, and 900 km at these three times. The values of the specific 
entropy in the convectively unstable gain region increase with time in all simulations. Model .v27/hcatl .00 exhibits a stagnant shock and only small deviations 
from sphericity. The average shock radius is secularly growing in model i'27/i,eatl .05 with slightly stronger neutrino heating and the shock is more aspherical. 
Model i27/[,eat 1 ■ 1 5 is on track to explosion and exhibits, at 154ms after bounce, a strongly deformed shock with a single large high-entropy bubble. 

(red). Red and yellow thus mark gas in the high-entropy gain 
layer, while cyan indicates the shock front and an iso-entropy 
surface at the edge of the protoneutron star While the shock 
appears nearly spherical in model .s'27/heatl OO, it is clearly 
deformed in model 527/heatl-05, and strongly so in models 
i27/heatl-10 and i27/heatl-15. One also notes that in the latter 
two models the highest-entropy gas is concentrated in the re- 
gion of greatest expansion while it is more evenly spread out 
in the other models. The shock deformation in these models 
is clearly dominated by low-£ modes, but there is still much 
smaller-scale structure in the form of protrusions caused by 
rising hot gas bubbles that push out the shock front at local 
scales. 



Figure [4]depicts colormaps of 2D x-z slices of the specific 
entropy in models .s27/heatl-00, s27/heatl-05, and .s27/heatl-15 
at 80, 115, and 154ms after bounce. Model .s27/heatl-10 is 
not shown, but is overall very similar to model i27/heatl-15. 



The evolution towards large shock radii, large-scale shock de- 
formation, and peak specific entropies of >20^b is obvious in 
the slices belonging to models 527/heatl-05 and .s27/heatl-15. 
In the latter, at 154ms, one notes a large high entropy area 
subtending an angle of ^30° and ranging from the gain ra- 
dius out to the shock, which has the overall greatest radii in 
this region. This is consistent with the volumetric view of this 
model at approximately the same time, shown in Fig.|3] 

In the bottom panel of Fig. [2] we plot the density-weighted 
average of the specific entropy m the gain layer (igain) (dashed 
lines; right ordinate). While the heating rates differ strongly 
between the models, their (.Sgain) remain very similar until 
~'100ms after bounce and (.Sgain) ^ 12^Bbaryon~'. The top 
row of Fig.|4]shows that the peak entropy reached in the gain 
layer is very comparable among the three displayed models at 
80ms after bounce. At 115ms and, in particular, at 154ms, 
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Figure 5. Top panel: Evolution of tiie mass accretion rate measured at tiie 
sliock (rigiit ordinate) and baryonic mass of the protoneutron star enclosed by 
the 10" 



'gem" 



^ density isosurface (left ordinate). Center panel: Evolution of 
the protoneutron stai' radius, defined as the location of the lO' ' gcm"^ point 
on the angle-averaged rest-mass density profile. The protoneutron star con- 
tracts as neutrino-cooling and deleptonizing material is settling on its surface. 
This is expected from ID and 2D neutrino radiation-hydrodynamics simula- 
tions (cf. ^B. Miiller et al._2012 bl, but is not captured by the simple MB08 
light-tjulb approach. Bottom panel: Evolution of the gray, angle-averaged 
neutrinosphere radii {R,j) as predicted by the leakage scheme. The well- 
known hierarchy R^^ > Rj^^ > R^^ is reproduced and the neutrinospheres fol- 
low the contraction of the proton eutron star as exp ected from full radiation- 
hydrodynamics simulations (e.g., |Janka et al.|2007} . 

the situation is different. The models with increasing shock 
radii and shock deformations develop large regions with spe- 
cific entropies in excess of 20^Bbaryon~' and large spatial 
variations. In model .s27/heatl-15, at 154ms after bounce, the 
expanding deformed shock has akeady swept up cold gas that 
now moves through the gain layer, leading to a decreasing 
(■Sgain) in this model. This is consistent with the decrease in 
(■Sgain) seen at the onset of explosion in the 3D and 2D simu- 
lations of Hanke et al. (2012 ). At the same postbounce time, 
the shock in model .s'27/heatl OO has receded to ^90 km and 
the distribution of specific entropy behind it is much more 
uniform than in the other models. Its average entropy con- 
tinues to increase despite the decrease in net heating (cf. top 
panel of Fig.|2|i. This is due to the combined effect of smaller 
shock radii and small deformation of the shock front. The av- 
erage specific entropy in model 527/heatl-05 also grows, since 
its net heating rate continues to stay high while its shock de- 
formation is still moderate and shock expansion has not yet 
become dynamical. 

3.2. Protoneutron Star, Neutrino Emission, and 
Thermodynamics of the Post shock Region 

The three-species leakage/heating scheme employed in our 
simulations goes beyond the MB08 light-bulb approach taken 
by many recent 3D hydrodynamic studies (e.g., ,Nordhaus| 
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Figure 6. Angle-averaged profiles of specific entropy s (top panel), tem- 
perature T (second panel), electron fraction (third panel) and the rest- 
mass density p (bottom panel) at representative postbounce times in model 
.s27/heatl 05. The data are taken from the AMR level encompassing the 
shock and, hence, do not extend to the full 180 km shown at early times. 
The smoothness of the curves is due entirely to the angle averaging. The 
profiles show the progressive deleptonization and contraction of the outer 
protoneutron star and the development of the high-entropy gai n layer as is 
expected from full radiation-hydrodynamics simulations (e.g., |Buras et al.| 
2006b a Lentz et al. 2012, B. Muller et al. 2012bi. Note, however, that our 
leakage/heating scheme tends to somewhat overestimate cooling and delep- 
tonization at optical depths of a few, leading to a dip in Yg and a local tem- 
perature minimum around 40 km. This temperature minimum is shown in a 
zoomed-in inset in the temperature panel. 



et al.| 2010HHanke et al.|[20T2l [Burrow s et al.|2012l |Murphy| 
et al. 2012; Dol ence et al.||2012[ ). These simulations use an- 
alytic cooling functions and neglect important protoneutron 
star cooling by . They also do not take i nto account changes 
of the electron fraction after bounce ( |Hanke et al.||2012] l 
or do so only via a parameterization of Yeip), which cannot 
account for the strong deleptonization in the region behind 
the shock due to electron capture on free protons. Neutrino 
heating is realized in these simulations by an analytic heating 
function with spatially and temporally constant neutrino tem- 
perature and luminosity. An important consequence of these 
approximations is that accreted material settling onto the pro- 
toneutron star cannot sufficiently cool, deleptonize and con- 
tract panke et al.||20T2t IBTMUller et a l. 2012a). This, in 
turn, results in too large shock radii and low advection speeds 
through the convectively unstable gain layer th at may artifi- 
cially favor the growth of convection over SASI ("Scheck et al. 
2008 i FogHzzo et al. 2006; B. Muller et al. 2012a). OurTeat 
age/heating scheme is designed specifically to overcome these 
limitations at little additional computational cost. We take 
into account cooling by i^^, v,,, and f^^, account for the change 
in electron fraction by v^, and emission and absorption. Our 
heating prescription uses the true and Vi, luminosities avail- 
able at a given position for heating (as computed by leak- 
age/heating at smaller radii) and the mean-squared neutrino 
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energies entering the heating rate are determined by assum- 
ing black body emission from the and neutrinospheres, 
taking the time-changing thermodynamic locations on these 
surfaces into account. 

W hile clearly not as sophisticated as recent gray multi-D 



(e.g. IScheck et al.'2008';'E. Muller et al.'2012';'Kuroda et al. 
2012|i or en ergy-dependent (e.g., Ott et al. 2008, M arek"^ 
Janka 2009; B. Mull er et al.| |2012b a, Taki waki et aLpUnj l 
neutrino radiation-hydrodynamics calculations, the goal of 
our approach is to capture the essential qualitative features 
correctly and reproduce quantitative results approximately. In 
the following, we investigate the extent to which our scheme 
lives up to its premise. 

In Fig. |5] we plot, for all four models, the time evolutions 
of the baryonic mass inside the 10"gcm~-' density isosur- 
face (top panel, left ordinate), the angle-averaged accretion 
rate measured outside the shock (top panel, right ordinate), 
the angle-averaged coordinate radius of the 10"gcm"^ den- 
sity isosurface (center panel), and the angle-averaged Ve, v^, 
and i^v neutrinosphere radii (where t^^ = 1; bottom panel). 
The evolutions of protoneutron star mass and radius, and of 
the accretion rate are very similar in all models. The ra- 
dius if the 10"gcm"^ isosurface, which we define as the 
s urface of the protoneutron star following |B. Muller et al.| 

!2012bl, shrinks from ^70 km early after bounce to 40 km at 
80 ms after bounce. At the same time, the enclosed baryonic 
mass increases from ^I.ISMq to 1. 55 Mq. If accretion sud- 
denly stopped completely at 180 ms, the gravita tional mass 
of the final, co ld neutron star would be ^ 1 AMq ( jLattimer & 
|Prakash|2001| l. The increase in mass and decrease in radius of 
the protoneutron star seen in our simul ations is qualitat ively 
consistent with the findings of B. Muller et al.|(|2012bl and 



|Buras et al. (2006a) for different progenitors. B. Muller et al. 
( 2012a[ l, who studied the s27 progenitor, do not show these 



quantities. Hence, a direct quantitative comparison is not pos 
sible. 

The angle-averaged neutrinosphere radii given in the lower 
panel of Fig. [5] show that the leakage scheme correctly re- 
produces the well known hierarchy R^^ > > R^^ of neu- 
trinosphere radii i n the postbounce preexplosion phase (e.g., 
[Janka et al.| |2007| l. One notes that models with larger /heat 
have slightly larger neutrinosphere radii. We attribute this to 
their somewhat hotter postshock regions, resulting in higher 
opacity. 

Model ,s-27/heatl05 is intermediate between model 
527/heatl-OO that fails to explode and models .s27/heatl-10 
and s27/heatl-15, which have rapidly increasing shock radii 
at the end of their simulations. We choose .s27/heatl-05 as 
our representative model and show, in Fig. [6] angle-averaged 
profiles of its specific entropy, temperature, electron fraction, 
and rest-mass density at 40, 80, 120, and 140 ms after bounce. 
The smoothness of the profiles is due entirely to angle 
averaging. The overall qualitative behavior of all quantities 
is as expected from more complex radiation-h ydrodynamics 
simulations (cf Buras et al. 2006b and Fig. 5 of |Dessart et al. 
[2006). The radial extent and specific entropy of the gain layer 
increase with time, while the changes in the specific entropy 
below ^40km simply reflect protoneutron star contraction. 
The latter is also well captured by the rising temperature 
at the protoneutron star edge, indicating compression. The 
strong deleptonization of the postshock region caused by the 
i/g neutronization burst shortly after bounce is still visible 
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Figure 7. Top panel: i/^ , Ue, and Uy luminosities as a function of postbounce 
time in models i27/i,e;,tl 05 (solid lines) and 4-27/[,eatl l5 (dashed lines) as 
representative examples of our model set. The Ue and Ue luminosities in 
model 4-27/heiitl l5 are somewhat smaller due to the strong charged-current 
absorption in this model. The inset plot shows the Ue deleptonization peak. 
Comparing the L^, shown here with those provided for the same progenitor in 
Fig. 8 of B. Muller et al. ( 2012a i demonstrates that our much more approxi- 
mate neutrino treatment still yields luminosities that agree within ~20% with 
the results of true radiation-hydrodynamics simulations. Bottom panel: Evo- 
lution of the mean neutrino energies {^i,.) in the same models, obtained via 
the assumption of black body emission at the respective neutrinospheres. Af- 
ter the transient very early postbounce phase, the usual hierarchy of mean 
neutrino energies is established. At >80— 100, ms after bounce, the evolution 
becomes qualitatively incoirect when (eiy,) an d {e,}-) surpass (ti/,). Compar- 
ison with the results o f|B. Muller et al.H2012a^ shows that this and the overall 
high predicted (e,,.) are an artifact of the leakage/heating scheme. 



in the profile at 40 ms after bounce. The outer postshock 
region re-leptonizes over time due to a slight dominance of i'^ 
over Dg absorption in the gain layer In the lower postshock 
region {R ^ 10-45 km), neutrino cooling and deleptonization 
continue and, as expected from more accurate neutrino 
transport calculations, a strong negative lepton gradient 
devel ops that m ay drive protoneutron star convection (e.g., 
[Dessart et al.|2006; Buras et al..2006a) . 

The top panel of Fig.[7]shows the total luminosities of Vg, v^, 
and Vx as predicted by our leakage/heating scheme for models 
527/heatl-05 and .s27/heatl-15, which we take as representa- 
tive examples. Differences between these models are minor 
and due to the greater heating in model 527/heatl-15. Since 
B. Muller et al. (2012a) provide these luminosities from their 

simulations in their Fig. 8, we can directly compare with 
flieir results. Their L^^ peaks at -385Bs"' (IB = 10^'erg), 
while ours peaks at ^365 Bs~' (a 5% dif ference). At 100ms 
after bounce, the [B. Muller et al^ ( [2012a| simulation suggests 

^62Bs~', while we find '--^68Bs~' in model .s27/heatl-15 
(~73Bs-' in model s27/heatl-05), a -10% (-20%) differ- 
ence. The and luminosities compare similarly well. 
The rather good agreement in total neutrino luminosities with 
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th e much more detailed ra diation-hydrodynamics simulation 
of B. Miiller et al. (2012a) suggests that our leakage/heating 
scheme captures the overall neutrino emission and its energet- 
ics in an acceptable way. 

The situation is different for the mean neutrino energies 
{en) shown in the lower panel of Fig. [t] We obtain esti- 
mates for {ey) of each species by assuming black body emis- 
sion from its neutrinosphere in the same way as for the mean- 
squared energies that enter the heating function (Eq.fTll. This 
kind of estimate is not reliable in the very early, highly dy- 
namical postbounce phase, but ^20 ms after bounce, the usual 
hierarchy of neutrino energies (e^^) > (epj > {e^^) is estab- 
lished analogously to the hierarchy of neutrino sphere radii 
(cf. lower panel of Fig. |5]l. This hierarchy is, however, bro- 
ken at times > 80- 100 ms, when the mean Ve and ener- 
gies exceed the mean energy of v^. This is clearly an artifact 
of our leakage/heating scheme and will not happ en in nature. 
It is also not found by B. Mii ller et al.| ( |2012a[ ). The reason 
for this incorrect behavior can be understood by considering 
the neutrinosphere radii plotted in the lower panel of Fig. 5 
and looking at the temperature and profiles shown in Fig. "E 
for model i27/heatl-05. At postbounce times > 80ms, one 
notices a global minimum in F^. around 40- 50 km. At the 
same location, a local temperature minimum develops. Both 
are related and caused by the inability of the leakage/heating 
scheme to establish a correct balance between emission and 
absorption at optical depths of a few. Unfortunately, the 
neutrinosphere recedes precisely into the local temperature 
minimum, while the Ve and v,, neutrinospheres sit in the lo- 
cal maximum at slightly greater radii. While the differences 
in temperature are not large, they are sufficient to explain the 
incorrect evolution of the {e^)- 

Comparing the values of {e^) predicted by t he leakage 
scheme with the results of |B. Miiller et al!] (|2012a) at times 
before the qualitative evolution becomes unrealiable, we find 
that the leakage scheme systematically overpredicts the mean 
energies. For example, at 50ms after bounce, in model 
i27/heatl.05, we find (e^^J -16MeV, (epj -18MeV, and 
(e^,) ~ 18.5 Me V . At t he same time the {e^) found by 
|B. Miiller et al.| ( |2012a| l ai-e, in the same order, 9.3 Me V, 
12.3 Me V, and 14MeV. 

In summary, the results shown in this section indicate that 
the leakage/heating scheme used in our simulations yields 
overally qualitatively correct thermodynamics/stratification in 
the postshock region and captures the integral neutrino emis- 
sion to within ^20% of fully self-consistent simulations. It 
fails, however, to yield reliable predictions for the mean neu- 
trino energies, in particular at later postbounce times. Since 
energy (and lepton number) absorption rates depend sensi- 
tively on neutrino energy, the leakage/heating scheme, at least 
in its present form, cannot be employed to make reliable pre- 
dictions of the spectrum of the emitted neutrinos or the com- 
position of explosion ejecta. 

3.3. SASI and Neutrino-Driven Convection 

The recent 2D radiation-hydrodynamics core collapse and 
postbounce sim ulations of the 527 progenitor carried out by 
B. MiiUer et al. ( j2012aj ) show a very clear and clean growth of 
a dominant periodic £ = 1 SASI mode. The relative amplitude 
(with respect to the average shock radius) of the t= \ mode 
saturates in their simulations at a very large ~45%, indicating 



a large-scale periodic dipole deformation. In their simulation, 
neutrino driven convection is only a secondary instability that 
develops in the non-linear phase, but may be connected with 
the saturation itself (Guilet et al. 2010). 

It is now interesting to ask if the SASI is the primary insta- 
bility driving asphericity in the ^27 progenitor also in our 3D 
simulations or if neutrino-driven convection dominates early 
on and possibly suppresses the growth of coherent SASI oscil- 
lations. It is furthermore interesting to study how the roles and 
prominence of SASI and convection depend on the strength of 
neutrino heating. It is evident from the discussion in §3.1 and 
Figs. |3]and|4]that deviations from sphericity develop at large 
scales m our models. We shall now take a more quantitative 
look at the development of this asphericity. 

3.3.1. Convection 

The local stability of a fluid element to convective o verturn 
is determined via the Ledoux criterion ( |Ledoux|1947| l, 



Cl = - 



dp\ 

dP) 



s.Yi 



dP 



p.Y, 



d.v 

dr J 



dP 

Wi 



dr 



(3) 

which, in the postbounce supernova case, takes into account 
radial gradients in specific entropy s and lepton fraction Y/ = 
Ye + Yu, ~Yc!^- For simplicity, we set 7/ := Y^, since our leakage 
scheme does not keep track of local neutrino fractions. This 
approximation may lead to quantitatively incorrect estimates 
of Cl in the protoneutron star where neutrinos are trapped or 
partially trapped. A fluid element is convectively unstable if 
Cl > 0. The linear growth time for convection from arbitrarily 
small perturbations is then given by the Brunt- Vaisala (BV) 
frequency. 



WBv = sgn(CL)' 



p dr 



(4) 



where we are following the definition of Buras et al. '(2006a'); 



Takiwaki et aL] ( 2012,) and where $ is the local gravitational 
potential and thus d'^/dr is the local gravitational accelera- 
tion. For simplicity, we approximate the gravitational accel- 
eration as -GM(r)r~^ assuming an angle-averaged spherical 
matter distribution in our postprocessing analysis. 



Foglizzo et al.| ( |200 6) pointed out that Eq. ^ is an insuf- 
ficient criterion for the development of large-scale convective 
instability in the postshock region. A small (linear) pertur- 
bation that could seed convection in the unstable gain layer 
is advected in towards the convectively stable cooling layer 
with the background flow. This advection may occur faster 
than the time it takes for convection to grow from the small 
perturbation. It is thus necessary to compare the advection 
timescale T^dv wit h the growth time for c onvection in the gain 
layer, Tconv ~ ^bv- Foglizzo et al. (2006 1 defined the quantity 



X = 



f^BV J Tadv 

--— t/r= 

"^conv 



(5) 



where is the radial velocity through the gain region. A 
small scale perturbation of magnitude Si„ entering the gain 
layer from above may at most grow by a factor exp(x) to 
^out = i^inexp(x ) during its advection through the gain layer 
( [Scheck et al.|p:008j . According to the linear analysis of 



13 



280 
240 
200 
160 
^ 120 



I T ] r I I I M r I I I I I 1 T f I I I 9 I T r E I 

: s27/heatl.05 COBv 

''shock,avg. 
- '^shockjtnin. 
^shcx:k,inax. 




20 40 60 80 100 120 140 160 
t - ^bounce [n^s] 



20 40 60 80 100 120 140 160 
i - ^bounce Ns] 



I 



1 T 



I 



-14.8 -7.4 0.0 0.77 



1.54 0,0 0.01 0.02 0.03 0.04 0.05 0.06 



440 

400E- 

360 E- 

320 
_280E- 
I 240 1- 
^ 200 r 

160 E- 

120 
80 
40 




I ' ' ' I ' ' ' I ' ' ' I ' ' ' I 

^— ^shock,inin. 
^— ^shock,ma>;. 




-14-9 



60 80 100 120 
f - ^bounce [^ts] 

-7.45 0,0 1.0 



20 40 60 80 100 120 140 
^ ~ ^bounce [ms] 



I 



2.0 0.0 0.01 0.02 0.03 0.04 0.05 0.06 



Figure 8. Coloraiaps showing the time evolution of the angle-averaged Brunt- Vaisala (BV) frequency wbv in units of ms"' (Eq.jJ] left panels) and anisotropic 
velocity Vjmiso in units of the speed of light c (Eq.l6] right panels) in models i27/heail-05 (top panels) and s27/hcatl l5 (bottom panels). Also indicated are the 
maximum shock radius (red curves), the average shock radius (blue curves), and the minimum shock radius (green curves). Note the different radial scales of the 
top and bottom panels. We mask out Va„iso and cjbv outside the average shock radius, where they are not reliable, since at most angles the radial region is actually 
outside of the shock. Shortly after bounce, the stalling shock leaves behind a negative entropy gradient, leading to cjbv > and thus convective instability, 
strongest at radii between 20 and 40km. Prompt convection develops quickly and is strong, as indicated by the large I'aniso in the light panels. Subsequently, 
convective instability and, as shown by the right panels, convection, develops in the gain layer and, after ~30-40ms or so, also at the edge of the protoneutron 
star core, due to the negative lepton gradient. Note that high I'aniso at late times prevails to significantly smaller radii than the inner radius of convective instability, 
indicating large asymmetries and undershooting of decelerating convective plumes. 
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Figure 9. Top panel: Foglizzo parameter x (Eq-[5j as a function of time 
after bounce. At times < 40 ms after bounce, tlie shock is still expanding 
and a quasi-stationary gain layer has not yet developed, x is not reliable 
in that phase. At later times, it stays consistently below the critical value 
of 3 suggested by Foglizzo et al. ( 2006 1 as being necessary for convection 
to develop from arbitrarily small perturbations. Note that stronger neutrino 
heating leads to greater since lubv is larger. Bottom panel: Density- 
weighted average of the anisotropic velocity (I'aniso) (Eq.pl in the gain layer 
inside the minimum shock radius. As in the case of x. uiis quantity is not 
reliable in the highly dynamic early postbounce phase. The early peak around 
10ms is related to prompt convection, which ebbs over ~30ms. Startin"at 
~40ms after bounce, when neutrino driving becomes efficient (cf. Fig.[2l, 
(i'aniso) increases nearly monotonically in a very similar way in all models. 
Only model i27/i,eatl 00, which has the weakest neutrino heating, deviates 
from this trend at times >100ms after bounce. 



[Foglizzo et al. ( 2006 1, X ^ 3 is required for convectio n to de- 
velop in the gain layer from small perturbations Si^. S check| 
|et al.| ( |20081 ) noted, then demonstrated, that the situation is dif- 
ferent if the seed perturbations 6in are sufficiently large so that 
the time integral of the buoyant acceleration becomes compa- 
rable to the advection velocity. In this case, the advected seed 
may grow into a buoyant plume and stay in the ga in layer in- 
stead of leaving it. The results of |Scheck et al.| ( [200 8 ) indicate 
that local seed perturbations of order 1%, e.g., in the upstream 
radial velocity, may already be sufficient to trigger convection 
even if x < 3. 

If convection does develop, a simple measure of its strength 
is the anisotropic velo city Vaniso, which we define, following 
ITakiwaki et aL] ( |20T2| , as 



(P) 



47r 



(6) 



where denotes an angle average at fixed radius. We com- 
pute Vaniso by introducing a spherical auxiliary grid onto which 
we interpolate the Cartesian coordinate velocity components 
and transform to obtain v^, vg, and v^p. We then integrate 
over Ait steradian at each radius r to obtain the various angle- 



averaged quantities. Vaniso is high in regions of large fluctua- 
tions in Vr and high non-radial velocities vq and v^. 

In Fig. [8] we present colormaps showing the time evolutions 
of radial profiles of the angle-averaged ljbv (left panels) and 
I'ani.so (I'ight panels) for models i27/heatl-05 (top panels) and 
i27/heatl-15 (bottom panels) as representative cases for mod- 
erate and strong neutrino heating. The qualitative evolution is 
the same in all models. 

Within milliseconds of bounce, a highly convectively unsta- 
ble region develops where the negative entropy gradient left 
behind by the stalling shock is strongest. As is evident from 
the v'aniso colormaps, a strong burst of prompt convection de- 
velops and smoothes out this entropy gradient within ^20 ms 
of bounce. The highly dynamical early phase of shock expan- 
sion and prompt convection is over by ^40 ms after bounce, 
when the shock has settled at ^100- 120km. At this time, 
the gain layer has developed and neutrino heating creates a 
negative entropy gradient and thus instability to convection 
between ^-^80 km and the shock. Also, deleptonization at the 
edge of the protoneutron star (cf. Fig. |6| creates a negative 
lepton gradient, driving protoneutron star convection, which 
sets in at 35 -40 ms after bounce and is clearly marked by a 
band of high Vaniso, spatially coinciding with the band of con- 
vective instability in the protoneutron star. 

Figure [91 s hows the evolution of the Foglizzo parame- 
ter X (Eq- [5j top panel) and the density-weighted average 
anisotropic velocity in the gain layer (bottom panel). At 
the early postbounce times at which prompt convection takes 
place, both quantities are poorly defined, since the shock ex- 
pansion is still rather dynamic and a quasi-stationary gain 
layer does not yet exist. This explains the large variations seen 
at early times in particular in x- Once the postbounce quasi- 
equilibrium in the postshock region is established, x settles at 
valu es between and 2 in all models, which is consistent with 
what B. Miiller et al. ( 2012a[ l found for the s27 progenitor. 

The X ?i 3 criterion proposed by 'Foglizzo et al. ( 2006 1 for 
the development of convection in the gain layer is never ful- 
filled in any of our models. Nevertheless, neutrino-driven 
convection does develop and becomes strong in all of our 
models. This is obvious from the radial Vaniso distribution 
shown in Fig. [8] As soon as the gain layer develops and cjgy 
becomes large, a broad region of high Vaniso appears and traces 
the region of instability. This is indeed neutrino-driven con- 
vection, as can be seen from the entropy slices in the top panel 
of Fig. [4] which show fully developed neutrino-driven convec- 
tion at ^80 ms after bounce. The development of neutrino- 
driven convection in the gain layer can also be inferred from 
the density-weighted average v'ani,so over the gain layer (bot- 
tom panel of Fig.|9]l. (va„iso) has an initial local maximum due 
to prompt convection and decreases as the latter ebbs only to 
increase again at > 40 ms after bounce when neutrino heat- 
ing in the gain layer becomes efficient and drives convection 
(cf.Fig.[2|). 



B. Miiller et al. ( 2012a l, in their axisymmetric simulation 
or the s27 progenitor, did not observe the development of 
neutrino-dri ven co nvection, in agreement with the prediction 



of Foglizzo et al. ([2006) that for x ^ 3 small perturbations 



are advected out ofthe gain layer before they can grow into 
buoyant plumes. While there are many technical dif ferences 
between the simulation of B. Miil ler et al.| ( |2012a| l and the 
ones presented here, the key difference relevant for the de- 
velopment of neutrino-driven convection in our simulations 
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is our choice of a Cartesian AMR grid as oppo sed to the 
spherical polar grid of the axisymmetric code of B. Miiller 



et al. ( 2012a| l. A spherical polar grid is ideal for tracking 
the spherically-symmetric collapse phase and the upstream 
flow outside the shock after bounce. Seed perturbations re- 
main minimal and neither prompt nor ne utrino-driven co nvec- 
tion grow in the simulation of B. Miiller et al. (2012a I. Our 
Cartesian AMR grid, on the other hand, leads to significant 
perturbations in multiple ways: (/) The Cartesian grid itself 
only imperfectly resolves spherical flow and perturbations of 
at most order dx/R^\^ac)L, where dx is one linear computational 
cell size, are generated locally at the shock front. (//) Also due 
to its rectangular nature, the grid has ^ = 4,m = 4 symmetry, 
which leads to buildup of numerical noise primarily in modes 
with £ = 4, m = {-4, 0,4}. (///) In our AMR setup, the shock is 
formed on the finest grid and then is allowed to pass through 
two mesh refinement boundaries before it reaches the grid that 
will track its subsequent evolution. The crossing of AMR 
boundaries causes large perturbations in the shock front that 
are also of ^ = 4 character (/v) The AMR grid that tracks the 
shock front expands whenever the shock expands. The AMR 
boundary must constantly be filled via interpolation from the 
next coarser grid, which also introduces noise. Points 
are true for any code using a Cartesian grid, e.g., the CASTRO 
code used in the rece n t simulations of [Burrows et al.| ( |2012[ ) ; 
Murphy et al] ([2012]); [Dolence et al] ( |2012[ ), whileOn)-(iv) 
are due to our particular approach in Zelmani, which may 
or may not be different from what is done in CASTRO and 
other codes. 

The £ = 4,m = -4,0,4 grid modes imprint themselves on all 
fluid variables, and we find that local variations in, for exam- 
ple, the radial velocity with respect to its angle-average are of 
order 1 - 10% behind the shock in the first few milliseconds 
after bounce. The shock shape itself remains relatively spher- 
ical, as can be seen in the bottom panel of Fig. [TO] which 
shows the root-square-sum A4 of the normalized T = 4,m = 
{-4, ... ,4} components of the shock front in all models (see 
Eqs.[8]and|9]l. A4 has its maximum of ^1.4% at ^lOms after 
bounce. 

The large-amplitude perturbations present in the early post- 
bounce flow are more than sufficient to overcome advection 
and seed prompt convection, which grows within milliseconds 
of bounce in our models (cf. Figs. [8]and|9]l. Neutrino-driven 
convection is, in turn, seeded by the turbulent flow of prompt 
convection and by additional, though much smaller magni- 
tude, noise coming from interpolation at the AMR boundary 
and from the Cartesian representation of the spherical shock. 

3.3.2. SASI 

Convective overturn, first prompt, then neutrino-driven, de- 
velops early on in our simulations and appears dominant. We 
can, however, not yet exclude growth of the SASI. The con- 
ditions for SASI growth are very different from those for con- 
vection. Any standing accretion shock is unstable to the SASI, 
with £ = l,m = Q,zLl m odes being the most unstable (e.g., 
Guilet & Foglizzo|2012| ). The linear growth rate of the SASI 
can be expressed as 



WSASI ■ 



InlQI 



(7) 



where Q is the cycle efficiency, defined as the amplification 
factor of perturbations in each advective-acoustic cycle, and 
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Figure 10. Time evolution of tlie normalized root-square-summed spherical 
harmonic mode amplitudes of the shock in all models for each £ in {1,2,3,4}: 
Al, A2, A3, and A4 (Eq,[9}. Note the dominance of the £ = 4 perturbations 
from the Cartesian grid at early times. The A[ amplitude becomes dom- 
inant ~ 40 — 50 ms after bounce and shows oscillatory features in models 
.v27 /heat 1 00 and i27/hciit 1 05. Models sllfi^ca 110 and .v27/heat 115 have no 
obvious oscillatory behavior of A 1 , but develop large non-oscillatory ampli- 
tudes at late times when the shock in these models reaches large radii. 



Tcyc is the duration of a cycle (see, e.g., Scheck et al. 2008 |f or a 



detailed discussion). Qualitatively, r^yc depends on the radius 
at which the shock stalls and on the timescale for advection of 
entropy/vorticity perturbations between shock and protoneu- 
tron star edge. A smaller shock radius and shorter advection 
time will thus lead to a smaller r^yc and faster SASI growth. 
Strong neutrino heati ng, as po inted out by Yamasaki & Ya-| 



mada (2007) and Scheck et al.| ( |2008 ), increases the buoyancy 
in the gain layer and leads to both larger Q and shock oscilla- 
tion frequencies (connected with Tcyc), while the growth rate 
is not strongly affected. 

Tell-tale signs of SASI growth are oscillatory deformations 
of the shock front. We look for evidence for the SASI in 
our simulations by decomposing the shock surface Rshockid, <P) 
into spherical harmonics: 



O-hn - 



(_l)l'"l 

^/47r(2£+l) 



R,bock{e,4>)Y^\9,<t))d^l. (8) 
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Figure 11. Normalized £= l,m = {—1,0, 1} mode amplitudes aim/'^oo of the shock front plotted on a linear scale (top panels) and their absolute values plotted 
on a logarithmic scale (bottom panels) for models .s27/heatl-00 (left panels) and s27fi^^^flA5 (right panels). Model .s27/heai 1 .00 shows a clear exponential growth 
of oscillatory modes, b ut satur ation occurs at amplitudes that are about an order of magnitude smaller than in the 2D simulation of the same progenitor cariied 
out by ,B. Miiller et ar|j2012a) . Model i27/iicat 115, which has strong neutrino heating and intense neutrino-driven convection also shows some oscillatory £ = I 
mode growth, though at a longer oscillation period, lower saturation amplitudes, and without a well defined exponential growth phase. 



Note that aoo corresponds to the average shock radius and that 
the definition of the aim used here gives individual ampli- 
tudes that are a factor of (2^-1-1) smaller than the definition 



for flfl) used by |B. Miiller et aL (2012a I in the axisymmetric 
case, but at each there are (2^-1- 1) more modes in our case. 
The Yj," are the standard real spherical harmonics (e.g., Boas 
2006 ), which we use with the normalization factors given m 
Burrows et al.| ( |2012j ). We also define the quantities as the 
root-square-sum of the a£,„ for a given £ normalized by the 
average shock radius ooo^ 



Af = 



\ 



7m 



(9) 



I Fig. [To] we present in four panels, from top to bottom, the 
; evolutions of the A1-A4 amplitudes of the shock front 



In] 

time 

in all four models. In the first ^20 ms after bounce, the ini- 
tial £ = 4 deformation due to our Cartesian grid imprints itself 
onto the shock front and the A4 amplitude is dominant. Subse- 
quently, the other modes grow. For SASI growth, the expec- 
tation is that the £= l,m = {-1,0, 1} modes have the fastest 
growth rate and have oscillatory behavior, which should be 
reflected in the Ai amplitude. In models 527/heatl-OO and 
i27/heatl-05, A I indeed is the fastest growing amplitude and 
shows the expected oscillatory behavior throughout the sim- 
ulated postbounce interval, suggesting the presence of the 
SASI. However, the maximum value of Ai reached is ^0.04, 
which is an order of magnitude smaller than what was re- 
ported by B. Mi iller et al.| ( |2012a| l for their 2D simulation of 
flie s21 progenitor. 

For the two models with stronger neutrino heating. 



527/heatl-lO and .s27/heatl-15, the situation is different. Their 
Ai and A2 amplitudes hover around very similar small 
values without obvious oscillatory behavior until ^ 100ms 
after bounce, when large-scale deviations from sphericity 
(cf. Fig. [3]l lead to strongly growing amplitudes in all £. 
Ai is the dominant amplitude and reaches ^0.1 in model 
i27/heat 1-10 and about 0.03 in model i27/heat 1 . 1 5 at the end of 
its simulation. It is interesting to note that model .s27/heatl -05, 
which has a positively trending shock radius at the end of its 
simulation, has clearly growing A2, A3, and A4 amplitudes at 
late times, while Ai remains the dominant mode with stable 
amplitudes near 0.03. 

For further insight into the nature of the observed mode evo- 
lution, we plot, in Fig. 1 1 the individual ^ = 1 , m = {-1 , 0, 1 } 
normalized mode amplitudes fli,„/aoo in linear (top panels) 
and logarithmic scale (bottom panels) for models .s27/heatl -00 
(left panels) and i27/heatl-15 (right panels). The former 
model has the weakest neutrino heating and least vigorous 
neutrino-driven convection of all our models while the lat- 
ter model has the strongest heating and most vigorous con- 
vection. All £ = I modes in model 527/heatl-OO show a clear 
oscillatory behavior and, importantly, an exponential growth 
phase between ~20 and ~^80ms after bounce can be made 
out. However, saturation occurs at low fli„,/floo ^ 0.01 for all 
modes. As noted before, this is an order of magnitude s maller 
than fou nd in the axisymmetric simulations of B. Miille r et al.| 
( 2012a) ) (in which neutrino-driven convection did not develop 
as a primary instability). 

Interestingly, some of the fli,„/floo modes in model 
■s27/heatl-15 do exhibit oscillatory behavior, though with 
larger periods than in model s21 /heat 1 00. This is expected for 
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SASI growth under the influence of strong neu trino heating 
( [Yamasaki & Yamada|2007| [Scheck et al.|2008| l. The growth 
also saturates more quickly at amplitudes that remain a factor 
of ^2 smaller than in model 527/heatl-OO until ^100 ms after 
bounce, when the mode growth becomes non-oscillatory. It is 
not possible to unambiguously and clearly identify a phase of 
exponential growth of the ai,„/floo modes in this model. 

In summary, there is clear evidence for SASI growth in our 
models. It is strongest in the model with the least neutrino 
heating and weakest neutrino-driven convection. It is weakest 
in the model with the most neutrino heating and the strongest 
neutrino-driven convection. However, even in the model in 
which SASI growth is strongest, the SASI saturates at ampli- 
tudes that are an order of magnitude s maller than in the 2D 
simulation of |B. Miiller et al.| ( |2012a| l, which did not have 
any neutrino-driven convection. These observations suggest 
that 3D neutrino-driven convection is indeed detrimental to 
the development of large-amplitude SASI. This confirms the 
findings of Scheck et al. (2008) and IGuilet et aL| (2010). Fur- 
thermore, our results show that both instabilities can coexist 
and grow at the same time, but even if convection is sup- 
pressed (a case we cannot study in our 3D Cartesian AMR 
code), the nearly equal splitting of the l= \ power across the 
three azimuthal m modes in 3D, will likely reduce the mag- 
nitude of deviations from sphericity that can be driven by the 
SASI alone. Moreover, the SASI, once it has reached its non- 
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hnear phase, wiU trigger neutrino-driven convection (|Scheck| Figure 12. Top panel: Ratio of advection and tieating timescales r,dv/T},eat 



let al.|2008||Guilet et al.|20T0l|B. MuUer et al.|2012a| l, which 
may very well become the dominant instability, in particular 
if neutrino heating is strong. 

3.4. Criteria for Neutrino-Driven Explosions 

The simulations presented here end before an explosion 
is fully developed in any of our models. Nevertheless, in- 
teresting trends can be observed. Models .s27/heatl-10 and 
527/heatl l5 have strongly positively trending shock radii at 
the end of their simulations. The shock in model ,s'27/heatl 05 
also expands at late times, but the development of an explo- 
sion is definitely more marginal. Model i27/heatl 00 has a 
receding shock and thus a rather negative prognosis regarding 
explosion. 

A variety of criteria for neutrino-driven explosions have 
been discussed in the literature and it is interesting to see how 
the trends observed in our models compare with what is ex- 
pected from theory and other simulation results. 

From the bottom panel of Fig. |2] we find that models with 
stronger neutrino heating and, thus, more vigorous neutrino- 
driven convection have systematically more mass in the gain 
layer (Mgain) that can absorb neutrino energy. The low- 
amplitude SASI seen in our models, which is strongest in 
models with weakest heating and convection, does not appear 
to have any positive effect on Mgain in our simulations. In 
models that are trending towards explosion, Mgain increases 
as shock expansion sets in. This is consistent with previous 
work giving the most optimistic prognosis for models with the 



as a function of time after bounce in our models (left ordinate) and -Epos, the 
volume integral over positive values of the total specific energy in the gain 
region (right ordinate). Tadv/'Tieat 1 is considered to be a condition for 
runaway explosion. It is satisfied by all of our models with optimistic out- 
look. Models i27/heai 1.15, i27/[,eat 1-10, and s27/i,eat 115 reach the threshold 
at ~100, ~115, and ~142ms after bounce, respectively. Roughly ~20ms 
later, these models are beginning to develop regions with positive total en- 
ergy, which may be interpreted as the onset of explosion. Bottom panel: 
Maximum of the ratio of the angle-averaged squared speed of sound to the 
angle-averaged squared escape velocity. According to the antesonic condi- 
tion of Pejcha & Thompson (2012), no solution for a spherical stationary 
accretion shock exists for max(c^/vgjj.) > 3/16 0.19 and an explosion is 
expected to set in in models that surpass this value. 

when the most optimistic models ac tuall y move to somewhat 
smaller i 



■J gain / 



(cf. the discussion in §3.1 1. Thus, in agreement 



with Hanke et al. (20121), the average entropy in the gain layer 



is not a good indicator for a model's potential for explosion. 

A criterion frequently used to diagnose neutrino-driven ex- 
plosions arises from the comparison of the timescale for neu- 
trino heating Theat and the advection times cale Tadv f or material 



to pass throught the gain layer ( Burr ows & Goshy 1993[ Janka| 
200T] [Thompson et"aL] [2005 ; Murp hy & Burrcws||2008p nf 



greatest Mgain (s-g-^ [M urphy & Burrows 2008 ; Scheck et al. 
|2008[[B7Muller et ai:|20I2bt|Hanke et al.|20i2| ). 

Also shown in the bottom panel of Fig. J2] is the density- 
weighted average of the specific entropy in the gain layer 
((■Seain))- All models, trending towards explosion or not, ex- 
hibit the same (igain) evolution until ~130ms after bounce. 



heating is faster than advection through the gain layer, then 
a fluid parcel entering the gain region may absorb sufficient 
energy to reach positive total specific energy and thus become 
unbound. For Tadv/Theat ^ 1. shock expansion should set in, 
further increasing Tadv and thus leading to positive feedback 
and runaway expansion. 

In our simplified analysis, we set Theat = |£^gain|/2net, where 
2net is the net integral heating rate in the gain layer and jiSgainl 
is the volume integral of the (Newtonian) total specific energy 
of material in the gain layer. There are a variety of possible 



definit ions for Tadv (cf. the discussions in Murphy & Burrows 
2008.;,Mare"k & Janka|20()9l|B. Muller eral..2012bi i. Here, we 
use the definition Tadv = ^/^gain> where Mgain is the mass in 
the gain region and M is the accretion rate through the shock. 
Note that this definition is different from what we use in the 
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computation of the Foglizzo x parameter (Eq.|5]l. 

In the top panel of Fig. [12] we plot Tadv/Theat as a func- 
tion of time after bounce forall of our models (left ordinate). 
The behavior is as expected: the two models .s27/heatl-15 
and 527/heatl lO, which are strongly trending towards explo- 
sion reach Tadv/Theat ^ 1 already at ~'100 and ^115 ms af- 
ter bounce. The marginal model i27/heatl-05 also shows in- 
creasing Tadv/Theat, which rcachcs 1 at ^ 142 ms after bounce. 
There is, however, no hope for model i27/heatl-00, where 
Tadv/Theat always remains below ^ 0.5. 



Also shown in the top panel of Fig. 12 is Epos (right ordi- 
nate), the integral energy of unbound material (with positive 
total specific energy). When Tadv/Theat ^ 1.4 in our models, 
material starts to become unbound and Epos grows rapidly. 
However, at the end of our simulations, it is still far away from 
the energy needed to unbind the entire envelope and lead to a 
canonical ^ IB core-collapse supernova explosion. 



Finally, in the bottom panel of Fig. 12 we plot the time 
evolution of the maximum of the ratio oTthe angle-averaged 
square of the speed of sound (c^) to the angle-averaged square 
of the escape velocity, which we approximate as {v^.^^.) w 
2GM{r)/r, where M(r) is the enclo sed baryonic mass. This 
ratio is interesting, since Pejcha & Thompson (2012) have re- 
cently derived the antesonic condition. 



max 



> 



3 

16 



:0.19 



(10) 



beyond which no solution for a stationary spherically sym- 
metric accretion shock exists, marking the transition to ex- 
plosion. While the expanding shocks in our models are far 
away from sphericity, we find values of max({c^) / (vl^^)) > 
0.2-0.22, which is consistent with the expectation of Pejcha 
& Thompson (20121. Model i27/heatl-00, which has the most 
pessimistic outlook, does not reach max{{cl) / (vl^^)) > 0.19, 
while the margi nal model, jf27/heatl-05 d oes. The prognosis, 
according to Pejcha & Thompson ( 2012| is thus similar to the 
one based on the Tadv/Theat > 1 runaway condition. 

3.5. Gravitational Wave Signals 

Besides the neutrino signals already discussed in ^3.2[ grav- 
itational waves (GWs) are the only other direct probe of the 
processes occurring in the postshock region and in the pro- 
toneutron star. The overall GW signature of core-col laps e su- 
pernovae has been reviewed in detail by Ott ( 2009 1 and Ko- 
[take (201 1 ) and we refer the interested reader to these reviews 
for an in-depth discussion of the various potential GW emis- 
sion processes and their underlying physics. 

GW observations of the next galactic core-collapse super- 
nova could provide important insight into the role and rel- 
evance of multi-dimensional fluid instabilities, rotation, the 



structure of the protoneutron star, and the nuclear EOS (|Dim- 


mebneier et al. 


2008 


IMarek et al.|2009||Yakunin et al.|2010[ 


Murphy et al. 


2009 


Rover et al. 2009 Ott 2009 Kotake 


201 l|l. Recently, Lo 


gue et al.|(|2012') carried out a proof-of- 



principle study, demonstrating that Bayesian inference allows 
to select between different explosion mechanisms for a galac- 
tic core-collapse supernova. The reliability of this depends 
on the availability of robust waveform predictions from sim- 
ulations. Most currently available core-collapse supernova 
waveforms come from 2D simulations (as summarized by 



Ott] 2009 Kotake et al. 201 li, which can predict only one 
of the two independent polarizations. In the context of nonro- 
tating or slowly rotating neutrino-driven core-collapse super- 
novae, only very few waveform predictio ns from 3D simula - 
tions without symmetry constraints exist, [^yer et al.| ( |2004[ ), 
carried out Newtonian 3D smoothed-particle hydrodynamics 
simulations with gray flux-limited diffusion neutrino trans- 
port and studied the GW emission from matter motions and 
asymmetric neutrino emission up to ^80ms after bounce in a 
variety of different precollapse configurations with an d with- 
out initial rotation and large-scale asphericities. Kotak eTt al.| 
(2009, 2011) performed Newtonian 3D hydrodynamic simu- 
lations with a light-bulb scheme (similar to MB08, but with 
a better approximation to changes in Yg). They used ana- 
lytic initial conditions, a fixed accretion rate and a fixed in- 
ner spherical boundary at 50 km, but were able to evolve for 
^500 ms and studied the GW emiss ion from matter dynamics 
and asymmetric neutrino emission. Scheidegger et al. ( 2010| l 
performed full 3D Cartesian (without inner boundary) New- 
tonian collapse and postbounce simulations of a slowly rotat- 
ing progenitor with neutrino leakage (but no heating). They 
employed a monopole approximation for gravity with rela- 
tivistic corrections and evolved to 100 ms after bounce. Re- 
cently, E. Miiller et al. (2012) presented Newtonian 3D post- 
bounce simulations with GR corrections to the monopole term 
of the Newtonian potential. They used a time-dependent inner 
boundary that contracts from 60-80kmto 1 5-25km over 1 s 
following the prescription of |Scheck et al.| ([2008), but were 
able to evolve multiple progenitor models for > 1.2 s using 
a ray-by-ray gray two-species approximate transport scheme 
(neglecting 1/^) and imposed neutrino luminosities at the inner 
boundary. They extracted and studied in detail the GW emis- 
sion due to matter dynamics and anisotropic neutrino emis- 
sion. 

While the simulations presented in this study do not have 
the more sophisti cated neutrino transport treatment of ]R] 
Miiller et al. ( |2012[ |, they do not have an artificial inner bound 



ary with imposed core neutrino luminosities, are carried out 
in full GR, and include the cooling due to h'^ emission from 
the protoneutron star It is, hence, worthwhile to study the 
GWs emitted by our models. We restrict ourselves to GWs 
from the dominant accelerated quadrupole matter motions and 
ignore GWs from asymmetric neutrino emission. The ratio- 
nale for the latter is that our simple leakage scheme is un- 
fit to give a reasonable estimate for the true neutrino radia- 
tion field anisotropy leading to GW emission. Moreover, as 
demonstrated by previous work ([Ko take et al.'2009' '201 1'; 'E] 
[Muller et al. 2012; Marek et al.p?J 09, Yakunin et al. 201017 
GW emission due to asymmetric neutrino emission occurs at 
too low frequencies to b e relevant for earthbound detectors 
such as Advanced LIGO ([Waldman (LIGO Scientific CoUab 



|oration)"2011'i 'Shoemaker"2010l, Advanced Virgo ([Accadia 



et al. (Virgo Collaboration) 201 1 ), and KAGRA ( |Somiya (for 
[lie KAGRA collaboration) |20 12) . 

We employ the quadrupole approximation for extracting 
G Ws from our sim ulations and use the expressions detailed 
in |Ott et aL] ( |2012) l. In principle, we could extract the grav- 
itational w aveforms directly from the spacetime, but the re- 
sults of Rei sswig et al.| ( [20lT| suggest that the quadrupole ap- 
proximation is very likely sufficiently accurate for stellar col- 
lapse spacetimes with a protoneutron star. The full observer- 
angle independent GW signals for all models are available 
for download from http : / /www . stellarcollapse^ 
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Figure 13. Left panel: Gravitational wave polarizations h+D and hxD (rescaled by distance D) of model i27/|,eat 1 05 as a function of postbounce time seen 
by and observer on the pole (6 = Q.ip = 0; top panel) and on the equator (8 = 7r/2,ip = 0; bottom panel). Right panel: The same for model s27/i,(;atl l5- Both 
models show a burst of gravitational waves associated with large-scale prompt convection developing shortly after bounce. Subsequently, gravitational wave 
emission comes from aspherical flow in the gain layer, in the outer protoneutron star, and from descending plumes of material that are decelerated at the edge of 
the protoneutron star. The gravitational wave signals are trending towards higher frequencies with time. 

convection in the two models are quite different, but the over- 
all amplitudes agree well, but peak in different viewing direc- 
tions. The subsequent evolution of the GW signals is similar 
in both models, both polarizations, and both observer posi- 
tions. After an intermittent quiescent phase, GW emission 
picks up again at times >80ms after bounce when aspherical 
dynamics becomes strong throughout the entire postshock re- 
gion (cf. Fig. [8]l. In this phase, the GW emission transitions 
to higher frequencies, indicating that emission from deceler- 
ation of downflows at the steep density gradie nt at the edge 
of the protoneutron star (as first pointed out by |Murphy et al.| 
|2009| l and convection in the protoneutron star play an mcreas- 
mg role. While both models have expanding shocks at the end 
of their simulations, the shock acceleration has not become 
sufficiently strong to lead to an offset in the GW signal (GW 
memory) seen in oth er work that followed exploding models 
to later times (e. g.. Murp hy et al.|20 09', 'Yak unin et al.|2010[ 
|E. Miilleret aLp 012, Kota ke et al.p 009 20lTf^ 
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Figure 14. Characteristic spectral strain spectra /ichar(/)/ of all four 
models at a distance of lOkpc compared with the design noise levels \/S( f ) of 
Advanced LIGO in the broadband zero-detuning high-power mode (aLIGO 
ZD-HP), KAGRA, and Advanced Virgo in wideband mode (AdV WB). 

[org/gwcatalog . 



The peak GW strain amplitudes reached in our models are 
from prompt c onvection and go up to ^20cm (^6.5 x 



10^ ^ at lOkpc). 



and 



Fryer et al. 



Scheidegger et al. (20T0| found \h\D 
(|2004 found 1^ 



^12cm, 



10cm 
but we note that 



In Fig. 13 we plot the h+ and hx polarizations of the GW 
signal (rescaled by distance D) for model 527/heatl.05 (left 
panel) and model .s27/heatl.l5 (right panel) as seen by ob- 
servers on the north pole (9 = Q,ip = Q; top panels) and on 
the equator (6 = Tr/2,(p = 0; bottom panels). The GW signals 
emitted by the other models are very similar and not shown. 
The early emission sets in ^ 10ms after bounce and is due to 
prompt convection that dominates the aspherical dynamics in 
the early postbounce phase, but has decayed by ^40 ms af- 
ter bounce. The GW signal from convection and other fluid 
instab ilities is of stochastic nature (cf K otake et al. 2009) |Ott| 
|2009| l and its time series cannot be predicted exactly. The GW 
signal of prompt convection, since it is emitted within mil- 
liseconds of bounce by the strongest first few overturn cycles, 
is particular sensitive to the perturbations seeding prompt con- 
vection. Note that the time series of h+ and hx from prompt 



the GW signal will depend on the strength of prompt convec- 
ti on, which is different fro m m odel to model. The a pproaches 
of |E. Muller etal] ^MT) and Kotake e t al.| ( |2009||2011| ) do 
not allow them to study prompt convection. The typical am- 
plitudes reached in the preexplosion phase are ~'3cm (^10~^^ 
at lO kpc). This is compara ble to, but somewhat larger than 
what E. Muller et al. (2012 1 found in the preexplosion phase 
of their models. This may be due the different progenitor 
models used and/or to the rather large inner boundary radius 
of their models in the preexplosion phase. Our typical \h\ are 
also quantitativel y consistent with the findin gs of the simpler 
3D simulatio ns of Scheidegger et al. ( 2010| l and Kotake et al. 



( 20091 2011 1, but are a fact or of a few smaller than predictions 
from 2D sirnulations (e.g. , |Marek et al. 12009) [Yakunin et al 
[20T0)[Mi^y et al.|2009| ). 



Figure 14 contrasts the angle-averaged characteristic GW 
strain spectra /ichai(/) (Flanagan & Hughes^ ^1998) of our 
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models with the broadband design noise levels of advanced- 
generation GW interferometers, assuming a source distance 
of lOkpc. The spectra are scaled with a factor of f~^l^ to 
allow one-to-one comparison with the detector one-sided am- 
plitude spectral noise density ^S{f), which has units of Hz''^^. 
Most of the detectable emission is within ^60- 1000 Hz and 
at essentially the same level of ^2-6 x 10"^"' Hz"''^^. A galac- 
tic event (at lOkpc) appears to be well detectable by the 
upcoming generation of detectors. All four models, while 
having distinct individual h+ and /ix time series that vary 
greatly in the time domain, exhibit essentially the same ro- 
bust spectral features, independent of /heat and the exact post- 
bounce time the individual models are evolved to. The low- 
frequency to intermediate-frequency emission is most likely 
due to prompt convection in the early postbounce phase, while 
the high-frequency peaks at ~'400Hz and ^900 Hz are most 
likely due to the deceleration of downflows at the protoneu- 
tron star surface and protoneutron star convection. A more 
detailed investigation of these features must be left to fu- 
ture work, since it would require multiple quadrupole inte- 
grals to isolate emission regions as done, e.g., by Murphy 
[et al. (2009). Our present simulations provide only one globaT 
quadrupole integral and we do not have sufficiently finely 
sampled output for a postprocessing analysis of the GW sig- 
nal. 

The total energy emitted in GWs is 3-4 x IO'^'MqC^ in 
all our models and about 50% of the emitted energy is due to 
the higher-frequency GW emission at later postbounce times. 
This finding is consistent with the 3D results of S cheidegger] 
l^al. (2010). E. Muller et al. (2012 ), on the other hand, found 
emitted GW energies of only ~'1O~"M0C^. Their models do 
not include prompt convection and emit most of their GW 
energy at frequencies below ~^400-600Hz. This, again, may 
be due to the different considered progenitor structures and/or 
to the inner boundary of their simulations. 

4. DISCUSSION AND CONCLUSIONS 

We have carried out four 3D general-relativistic core 
collapse and postbounce simulations of the 21 -M(^ solar- 
metallicity progenitor of [Woosley et al. (2002i, systemati- 
cally varying the rate of neutrmo energy deposition to study 
the effect of variations in neutrino heating on the 3D post- 
bounce evolution in general and on the standing accretion 
shock instability (SASI) in particular. These simulations nei- 
ther employed an artificial inner boundary nor did they make 
any symmetry assumptions or approximations for the gravita- 
tional field. The resolution of our simulations is nearly twice 
as high and we carried them out for nearly twice as l ong as 
the only previous 3D GR study of Kuroda et al. p012| l. 



For neutrinos, we used an energy-averaged (gray) three- 
species neutrino leakage/heating scheme in the postbounce 
phase, whose only free parameter is a scaling factor in the 
energy deposition rate. The leakage scheme captures the es- 
sential aspects of neutrino cooling, lepton number exchange, 
and neutrino heating as predicted by fully self-consistent ID 
and 2D neutrino radiation-hydrodynamics simulations. Im- 
portantly, our simulations do not suffer from the limitations of 
simpler analytic "light-bulb" heating/cooling schemes, which 
cannot capture the contraction and deleptonization of the pro- 
toneutron star and result in artificially large shock radii and 
overestimated advection times through the postshock region. 



The light-bulb approach, due to its simplicity and low com 
putational cost, is being employed in many contemporary 3D 
simulations (e.g., [ No rdhaus et al."2010', 'Hanke et al."2012| 
[Burrows et al.|2012t [Murphy et al. 2012 , Dolence et al. 20121 
Howe ver, as pointed out by Ha nke et al.| ( |2012 ) and B. Muller] 



et al. (2012a), light-bulb calculations may yield qualitatively 



mcorrect results for the postbounce hydrodynamics and the 
respective roles and relevance of neutrino-driven convection 
and the SASI. 

Our approach was designed specifically to avoid the prob- 
lems of the light-bulb scheme and provide a realistic post- 
bounce setting for more robust conclusions on the postbounce 
evolution and the role of hydrodynamic instabilities. At the 
same time, our leakage/scheme is still computationally much 
cheaper and simpler than the approximate gr ay or energy- 
dependent 3D neutrino transport schemes of IKuroda et aLI 
([20T2|), ITakiwaki et al.|([20 l2), and E. Muller et al. (2012f; 
Wongwathanarat et al. ( 2010) ). This affords us with the ability 
to carry out parameter studies with high numerical resolution 
as presented in this work for the 21 -Mq progenitor 



[B. Muller et al.| ( [2012a| l previously carried out an axisym- 
metric (2D) simulation of the same 21-Mq progenitor with 
their 2D GR radiation-hydrodynamics code. They found 
neutrino-driven convection to be suppressed due to the high 
postbounce accretion rate and, thus, short advection time 
through the convectively unstable gain layer The SASI is the 
primary instability in their simulation and seeds convection, 
which grows only as a secondary instability once the SASI 
has reached non-linear amplitudes. 

Our models show instead early and strong growth of con- 
vective instability. It is initially prompt, driven by the nega- 
tive entropy gradient left behind by the stalling shock. Sub- 
sequently, convection is driven by neutrino energy deposi- 
tion in the gain layer Neutrino-driven convection first mani- 
fests itself in small-scale local rising hotter and sinking cooler 
blobs of postshock material. In models with strong neutrino 
heating that are trending towards explosion, the small scale 
blobs combine over time to a few large, near volume-filling 
high-entropy regions whose expansion pushes out the shock. 
These large blobs lead to a low-^-mode dominated structure 
of the expanding shock. The shock, however, has a com- 
plicated substructure of protruding bumps caused by smaller- 
scale plumes that perturb it locally. Models whose shock ex- 
pansion becomes dyna mical surpass the runaw ay explosion 

criterion Tadv/rheat > 

and satisfy the antesonic 

(2012). Both criteria for explosion yield predictions consis- 
tent with the trends in our models. Interestingly, shortly after 
the Tadv/Theat ^ 1 Condition is met by our models, individual 
fluid cells behind the shock reach positive total energy, indi- 
cating the transition to explosion. 

While neutrino-driven convection is the fastest growing and 
dominant instability, our analysis suggests that all of our mod- 
els exhibit some growth of clearly periodic low-£ deforma- 
tions of the shock front. These are the tell-tale signs of the 
SASI, and, as expected from linear perturbation analysis, we 
find that the £ = l,m = {-1,0, 1} modes exhibit the fastest 
growth. However, our results also show that the saturation 
amplitudes of the oscillatory ^=l,m = {-l,0, 1} modes are, in 



la mical surpass the runaw' ay explosion 
1 ( [Burrows & Goshy|1993[Jja nka 2001) 
onic condition of 'Pejcha ^^Thompsdn] 



the best case, an order of magnitude smaller than in B. Miiller 
[et al.| ( [2012a[ l. The SASI remains a sub-dominant instab 
ity in all of our models. Furthermore, we find the SASI to 
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be strongest in the model with the least neutrino heating and 
the weakest neutrino-driven convection. Models with stronger 
heating and more vigorous convection have lower saturation 
amplitudes of the oscillatory modes, but develop large non- 
oscillatory deformations of £ = 1 , 2, 3 character that are caused 
by low-mode neutrino-driven convection and are unrelated to 
the SASl. 

Our simulations satisfy all the requirements laid out by [BT] 
|Muller et al.| ( |2012a| for the development of strong SASI m 
the 21-Mq progenitor: GR gravity, an EOS that results in 
a fairly compact protoneutron star, and the inclusion of all 
neutrino species and deleptonization of the protoneutron star. 
Yet, our resu l ts turn out to be very different from what |B7| 
IMtiller et al.| ( |2012a| l found. What is the root cause of this 
discrepancy? On the one hand, our simulations are 3D, split- 
ting, on average, the 2D £ = 1 SASl power across three az- 
imuthal m modes. This may explain lower saturation ampli- 
tudes, but cannot explain the early growth of neutr ino-driven 
convection th at is absent from the 2D simulation of IB. Miillerl 
|et al.| ( |2012a| l. On the other hand, - and, as we are convinced, 
more importantly - our simulations used a central Cartesian 
adaptive-mesh refinement (AMR) grid, which imparts pertur- 
bations of order of 1 - 10% onto the very early postbounce 
flow, seeding prompt convection. This, in turn, acts as seed for 
neutrino-driven convection in our models. The seed perturba- 
tions are sufficiently large for convection to develop despite 
the high accretion rate and correspondingly short advection 
time through the gain layer. Neutrino-driven convection be- 
comes dominant and limits the growth of the SA Sl, in agree- 
ment with the 2D work of Scheck et al. ( 2008| l. We expect 
any 3D simulation relying on 3D Cartesian AMR with sim- 
ilar resolution to have similarly large seed perturbat ions for 
convection. The recent 3D light-bulb si mulations of|Burrows| 
[etaL] ( [20T2l ); [Murphy et"al. ( .2012j ; ,Dolence et al.| ( |2012[ ) are 
all subject to these perturbations. 

The questi on of the magnitude of seed perturbations was 
not raised by 'B. Miiller et al. ( 2012a i, who used a spherical- 
polar grid that leads to only minute perturbations from the 
growth of numerical noise during co llapse. Is the almost per- 
fectly spherical postbounce state of B. Miiller et al. ( 2012a| l 
representative of nature or should one expect significant as- 
phericities to be present in the outer core? Some guidance 
on the size of perturbations induced by turbulent convection 
during late time burning in core-collapse supernova progen- 
itors is already available from the 2D and 3 D simulations of 
Meakin, Arnett, and collaborato rs (Meakin 2006 [ [Meakin &| 
|Arnett|2006||2007b[|Arnett & Meakin,2011) ). 

There are two important results from these multi- 
dimensional stellar evolution calculations that pertain to 
the expected density perturbation amplitudes in precollapse 
cores. First, 2D and 3D simulations of the oxygen shell 
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burning dominated phase in a 23-M0Star (Meakin & Arnett 



2007b]) have clarified the basic mechanism responsible for the 



origin of the fluctuations. In short, Meakin & Arnett (2007a) 
found that the root-mean-square (rms) density fluctuations are 
largest at the convective boundaries. By interpreting the dy- 
namics of the convective boundary layer in terms of g-modes 
excited by the turbulent convection, it was shown that the 
rms density fluctuation amplitude can be related directly to 
the background stellar structure and the Mach number of the 
convective flow, with 



where Mc is the rms Mach number of the convective flow, wbv 
is the Brunt- Vaisala frequency in the stable layer adjacent to 
the convection zone, v, is the sound speed of the gas, and g is 
the gravitational acceleration. The first term on the right hand 
side, which is very small, is relevant to the interior of the con- 
vection zone, where density fluctuations arise solely from the 
presence of velocity fluctuations in a nearly adiabatic layer. 
The second term is significantly larger and applies to the sta- 
ble layers bounding the convection zone, reflecting the exci- 
tation of fluid motions in these regions in the form of internal 
waves (predominantly g-modes). 



The 27 -Mq progenitor of Woosley et al. "2002' has a turbu- 
lent Mach number of ^0. 1 to 0.2 in the silicon burning con- 
vective shell overlying the core, and two peaks in wbv of im- 
portance: the peak corresponding to the inner edge of the ac- 
tive silicon burning shell (corresponding to the outer edge of 
the iron core), and a peak deeper in associated with the out- 
ermost extent of the now extinguished silicon burning core. 
Both peaks have values of v^^Bv/g of ^1, indicating that rms 
density fluctuations at these locations will be of order the tur- 
bulence Mach number of the convection, or ^ 10-20%. The 
spike in ojbv associated with the outer extent of the silicon 
core burning epoch will be accreted into the shock within 
^15ms of bounce, while the edge of the iron core will be 
accreted a little later, at ^60 ms after bounce. 

The second result from the multi-D stellar convection sim- 
ulations of Meakin and Arnett involves the interaction of nu- 
clear burning shells at late times. While the results on bound- 
ary layer fluctuations described above are considered to be 
robust by those authors, the presence of two or more convec- 
tive shells in close proximity, as found in late burnings stages, 
has been found to drive additional motion at the convective 
boundaries and correspondingly larger density fluctuation am- 
plitudes. In the most relevant case of a silicon burning shell 
around an iron core, the interaction between the silicon, oxy- 
gen, neon, and carbon shells were found to produce a dramatic 
increase in boundary layer distortion, eventually leading to a 
complete disruption and mixing of the multi-shell burning re- 
gion ([Meakin 2006; Arnett & Meakin 2011). This result is 
likely to be due, at least in part, to the inconsistency between 
the initial stellar model used (based on mixing length theory) 
and a more realistic turbulent convection as represented by the 
numerical simulation. Judging the robustness of these shell- 
interaction results, however, awaits 3D simulations since all 
of the multi-shell calculations performed to date have been 
restricted to 2D geometry which is known to result in exag- 
gerated velocities in regions of thermal convection. From this 
body of work, it would appear that the presence of density 
fluctuations with amplitudes of at least 1%, and possibly as 
large as 10 to 20%, should be expected in the material accret- 
ing into the shock at early postbounce times in a collapsing 
iron core. 

The fast growth of neutrino-driven convection in our cur- 
rent models is almost certainly caused by the large seed per- 
turbations from our Cartesian AMR grid. In 2D simulations, 
the growth of neutrino-driven convection may go along with 
SASl growth or, if not genuine SASl, then at least large-scale 
oscillatory low-£ deformations of the shock front {B. Mtiller| 
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et al. 2012a Burrows et al. 2012 Fernandez & Thompson 
2009a[ IScheck et al.||2008[ ). Our 3D models do not exhibit 
any large-scale oscillatory features. Rather, models evolv- 
ing towards an explosion develop non-oscillatory large-scale 
asphericities at late times and produce a globally aspherical 
explosion morphology without a need for SASl-driven £ = I 
deformations. This qualitative finding is in agreement with 
the results of the co nvection-dominated 3 D N ewtonian light- 
bulb calculations of |Burrows et al.| ( |2012| l and Dolence et al.| 
(2012^ . The late-time development of SASI-like oscillatory 
behavior seen in 2D simulations that are initially convection 
dominated (e.g., Marek & Janka 2009; B. Muller et al. 2012b]l 
may thus be an artifact of axisymmetry, but further work is re- 
quired to solidify this conclusion. 

The next galactic core-collapse supernova will reveal its in- 
ner workings by means of its neutrino and gravitational-wave 
(GW) signals. Both will provide key insight into the ther- 
modynamics and multi-D dyna mics of the protoneutron star 
and th e postshock region (e.g., |Ott||2009| |Lund et al.||2010| 
|2012| [O'Connor & Ott|2012| l. Whileour neuti'ino treatment 
is too simplistic to yield quantitatively interesting predictions 
of the neutrino signal, we are in a good position to study the 
GW emission from accelerated quadrupole mass motions in 
our models: For the first time, we have extracted GWs from 
full 3D GR collapse and postbounce core-collapse supernova 
simulations. We find a strong burst of GWs associated with 
early-postbounce prompt convection with frequencies around 
~100-200Hz, a subsequent almost quiescent phase, fol- 
lowed by higher-frequency (400- 1000 Hz) emission, whose 
amplitudes are dominated by the deceleration of undershoot- 
ing convective plumes at the edge of the protoneutron star If 
convection (prompt and/or neutrino-driven) does not develop 
early, the GW signal would not have a strong initial burst, but 
rather a slow rise to smaller amplitudes at later times, when 
the SASI becomes strong. This is a key difference and may 
allow GW data analysts to distinguish between convection- 
dominated and SASl-dominated postbounce evolution in the 
next galactic core-collaspe supernova. The design sensitivi- 
ties of advanced-generation GW detectors such as Advanced 
LIGO, Advanced Virgo, or KAGRA are likely to be suffi- 
cient to detect the collapse and neutrino-driven explosion in 
our 27-Mq progenitor throughout the Milky Way. While dif- 
ferent in detail, our results for the GW signature are gener- 
ally consistent with what was found for other progenitors in 



the 2D first-principles simulations of |Marek et al.J (i2009) and 

she 



|Yakunin et al.| ( |2010| l. Our GW signals have higher ampli- 
tudes and characteristic frequencies than predicted by the 3D 
simulations of E. Muller et al. (2012), who employed an artifi- 
cal inner boundary that was moved in according to an analytic 
prescription. 

There are a number of shortcomings and limitations of the 
simulations presented here that must be mentioned and can 
be removed only by future work. As is well kno wn and has 
been pointed out recently by Hanke et al. (20121 in the core 
collapse context, in 3D, turbulent power cascades to small 
scales. Low resolution in 3D may artificially keep power at 
large scales and may thus lead to an overestimate of the pos- 
itive effect of neutrino-driven convection. While our effec- 
tive angular and radial resolution in the postshock gai n layer 
is comparable to the highest resolution considered by jHanke 



et al. (2012 1, we agree with their assessment that understand- 
ing the resolution dependence of 3D results is of great impor- 
tance. We will carry out a resolution study in future work. 



The second major limitation of our simulations is our Carte- 
sian AMR grid and the fact that we must let the nascent su- 
pernova shock pass two mesh refinement boundaries before 
tracking its further evolution by AMR. This induces large 
perturbations leading to the growth of prompt and neutrino- 
driven convection in all of our models. These large and essen- 
tially unavoidable seed perturbations for prompt and neutrino- 
driven convection make it difficult to draw conclusions on 
which hydrodynamic instability dominates in the early post- 
bounce phase. This limitation is shared by other Cartesian 
AMR schemes and could possibly be avoided in future work 
by extending our spherical-polar grid blocks all the way into 
the protoneutron star core and using a single high-resolution 
Cartesian mesh only in the innermost few kilometers. 

A third major limitation of our work is the reliance on our 
simple gray heating/leakage scheme. While superior to the 
light-bulb approach, it cannot replace the energy-dependent 
neutrino radiation-hydrodynamics treatment that has proven 
to be crucial for reliable conclusions on the neutrino mecha- 
nism (e.g., |B. Muller et al. 2012b and references therein). The 
set of 3D general-relativistic hydrodynamics simulations pre- 
sented here required about ^20 million CPU hours to com- 
plete. Adding energy-dependent 3D neutrino transport will 
increase the computational complexity by an order of mag- 
nitude. Novel, highly efficient and scalable approaches to 
3D neutrino transport will be needed to address this problem 
(Sumiyoshi & Yamada 20 12[|Abdikamalov et al.|2012t[aiangl 
let al.,2012|,Radice et al.,2012[ ). 
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